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1.  Introduction 


This  report  is  concerned  with  the  upward  continuation  of  gravity  anomaly 
data  given  on  the  surface  of  the  earth.  A  number  of  computational  procedures 
to  be  presented  here  are  offshoots  from  the  studies  of  Cruz  (1985)  dealing  with 
modeling  the  external  gravity  disturbance  vector  field. 

Various  modern  possibilities  exist  to  model  the  external  gravity  field  of  the 
earth  for  the  realistic  case  of  free-air  anomalies  being  given  on  the  surface  of 
the  earth’s  visible  topography.  These  possibilities  fall  under  two  general  types 
of  modeling  approaches:  the  continuous  approach,  and  the  discrete  approach. 
In  the  continuous  approach  the  free-air  anomalies  are  assumed  to  be  known  at 
every  point  of  the  earth’s  surface,  and  the  Molodensky  problem  is  being 
solved.  In  the  discrete  approach  the  anomalies  are  known  only  at  discrete 
points  of  the  earth’s  surface,  and  the  Bjerhammar  problem  is  then  being 
solved.  In  this  report  we  have  used  the  discrete  approach,  specifically  least 
squares  collocation,  only  in  the  first  stage  of  data  processing,  for  the  purpose 
of  generating  an  optimal  set  of  mean  surface  free-air  anomalies  from  the 
originally  given  irregular  and  discrete  distribution  of  point  anomaly  data. 
After  this,  using  the  optimal  set  of  surface  data  and  concepts  from  the  con¬ 
tinuous  approach  we  generated  our  quantity  of  interest  which  is  the  upward 
continued  anomaly. 

The  simplest  conceptualization  of  a  solution  to  the  (continuous) 
Molodensky’s  problem  is  by  means  of  analytic  continuation  advocated  in  Moritz 
(1969).  The  external  gravity  anomaly  field  is  analytically  continued  to  a  level 
surface  which  may  be  entirely  above,  partly  above  and  partly  below,  or 
entirely  below  the  earth’s  surface.  Once  the  level  surface  anomalies  are 
known,  then  under  a  spherical  approximation  the  external  gravity  field  can  be 
generated  from  these  anomalies  using  classical  proc'  dures  for  data  on  a 
sphere.  A  general  procedure  for  analytical  continuation  is  by  means  of  Taylor 
series: 


Agp*  =  Ags  +  (Hp- 


Hs) 


+  2  ^(HP-  Hs)Z  " 


(1.1) 


where 

4t?s 

agp* 

Hs 

Hp 

<?Ag 

dH 


surface  free-air  anomaly,  defined  more  precisely  in  section  3 

anomaly  in  the  same  plumb  line  as  ags,  but  located  on  1  he 
level  surface,  to  which  the  surface  data  are  being  reduced 

elevation  of  the  surface  point  to  which  Ags  applies 

elevation  of  the  level  surface  of  4gp* 

vertical  gradient  of  the  gravity  anomaly  field 


If  the  level  surface  to  which  the  data  are  reduced  is  entirely  below  the  earth’s 
surface  the  analytical  (downward)  continuation  may  also  be  done  by  an 
inversion,  usually  by  successive  approximations,  of  the  classical  Poisson 


The  use  of  (1.1)  presents  practical  difficulties  because  the  computation  of 
the  required  vertical  gradients  of  the  gravity  anomaly  field,  even  just  the  first 
gradient  and  the  more  so  the  higher  order  gradients,  poses  rather  severe 
requirements  on  the  density  and  accuracy  of  gravity  data  (for  a  study  of 
numerical  evaluation  of  the  gradients,  see  Noe,  1980).  Therefore,  the  tech 
niques  that  we  have  used  in  this  report  have  a  common  motivation,  namely  to 
avoid  altogether  the  use  of  correction  terms  in  (1.1)  and  therefore  use  only  the 
first  term.  The  different  manners  in  which  the  correction  terms  are  avoided 
give  rise  to  three  methods  for  upward  continuation  which  we  numerically  tested 
and  compared  using  real  gravity  and  elevation  data.  The  first  method,  and  the 
crudest,  is  to  simply  drop  the  correction  terms  and  take  Agp*  to  be  equal  to 
Ags;  we  therefore  simply  insert  Ags  directly  into  the  classical  Poisson  upward 
continuation  integral  for  data  on  a  sphere.  This  procedure  will  obviously  be  in 
error  especially  in  areas  with  rough  anomaly  field  and  our  numerical  study  will 
provide  a  feeling  for  the  magnitude  of  this  error.  The  second  method  is  to 
drop  all  correction  terms,  only  after  the  terrain  correction  has  been  applied  to 
A  g  s .  The  application  of  the  terrain  correction  is  viewed  as  a  first  order 
attempt  to  reduce  the  surface  data  Ags  to  a  level  surface;  the  reduced  data  are 
then  inserted  into  the  Poisson  upward  continuation  integral.  The  third  method 
is  to  drop  the  correction  terms,  only  after  smoothing  the  anomaly  field  by  the 
subtraction  of  the  gravitational  effects  of  certain  shallow  topographic  masses  of 
assumed  density.  The  total  upward  continued  anomaly  field  is  then  the  sum  of 
two  fields:  one  generated  by  classical  Poisson  integration  from  the  residual 
anomalies  left  after  the  removal  of  topographic  masses,  and  the  other  is  the 

field  generated  by  integration  of  the  gravitational  effects  of  the  removed 

topographic  masses  themselves.  This  third  method  is  in  the  spirit  of  the 

•ni<  >v.  •  v ! :  •  •<  I  hi  t  !n  •  \  i  iid  .  f  "  '  and 

Forsberg  (see  Tscherning,  1979;  Tscherning  and  Forsberg,  1983;  and  Forsberg, 
1981,  p.  11). 

Other  concerns  of  this  report  include  the  use  of  spherical  harmonics  in 
anomaly  field  modeling,  the  use  of  Fourier  series  for  upward  continuation,  and 
the  application  of  studied  modeling  techniques  to  the  balloon-borne  gravity 
project  being  coordinated  by  the  Air  Force  Geophysics  l  aboratory  (Lazarewic/., 
et  ah,  1983). 


2.  Upward  Continuation  Formulas 


2.1  Spherical-Earth  Poisson  Integral 

Let  us  use  the  spherical  coordinates  r  (geocentric  radius),  ^  (geocentric 
latitude),  and  X  (geocentric  longitude).  The  anomalous  potential  T(r,  ?,  X) 
being  a  harmonic  function  in  space  has  surface  spherical  harmonics  Tn(R,  T,  X) 
attenuating  with  r~(n+1)  (Heiskanen  and  Moritz,  1967,  p.  35): 


Tn(r,  ?,  X)  =  (P)n+1  Tn(R,  *,  X) 


(2.1.1) 


The  surface  harmonics  of  the  gravity  anomaly  Ag(r,  ?,  X)  and  anomalous 
potential  T(r,  X)  are  fundamentally  related,  in  a  spherical  approximation,  as 
follows  (ibid.,  pp.  88-89): 


4gn( r*  x) 


X) 


(2.  1.2) 


The  last  equation  becomes  for  r-R: 


4gn(R.  ▼  ,  x)  =  Tn(R,  ?,  X)  (2.1.3) 

Substituting  (2.1.1)  into  (2.1.2): 

Agn-'r,  *,  X)  =  ^l-  (P)n+2  Tn(R,  ?,  X)  (2.1.4) 

Substituting  (2.1.3)  into  (2.1.4): 

4gn(r.  *,  X)  =  (p)n"2  Agn(R,  ?,  X),  (2.1.5) 

that  is,  the  surface  harmonics  of  Ag  attenuate  as  r-(n+2)i  The  upward 
continued  anomaly  Ag(r,  ?,  X)  is  found  by  summing  the  terms  in  (2.1.5).  After 
omitting  the  zero  and  first  degree  harmonics  as  customary: 

<o 

Ag:r,  ¥,  X)  -  Y^_  (?)n+2  4gn(R.  x)  (2.1.6) 

n "  2 


3 


The  space  domain  equivalent  of  (2.1.6)  is: 


Ag(r,  ▼  ,  X)  =  ^  |  J  K(t,  f)  Ag(R,  V,  X’)  dc7 

O' 

where 


t 


(2.1.7) 


K(t,  f) 


CD 


E 


(2n+l)  tn+2  Pn(cos  i>) 


(2.1.8) 


cos  =  sin  ?  sin  3>’  +  cos  <F  cos  <F’  cos  (X’  -  X) 


(2.1.9) 


A  closed  form  for  K(t,  V')  can  be  obtained  using  the  following  relation  (ibid.,  p. 
35): 


=  (2n+l)  tn+l  Pn(cos  i>)  (2.1.10) 

n-0 

where  D  (1  *-  t2  -  2t  cos  V')^  (2.1.11) 

Multiplying  (2.1.10)  by  t,  removing  the  zero  and  first  degree  harmonics,  and 
combining  with  (2.1.8)  we  get  the  closed  form: 


t2( 1  — t 2 ) 

K(t,  f)  =  — L  -  t2  -  3t3  cos  V  (2.1.12) 


Equation  (2.1.7)  with  (2.1.12)  is  the  same  as  equation  (2.160)  in  ibid.,  p.  90, 
and  is  in  fact  the  well-known  Poisson  integral  formula  for  the  space  domain 
upward  continuation  of  gravity  anomalies  given  on  the  surface  of  a  geocentric 
sphere.  Note  that  the  second  and  third  terms  of  the  integral  kernel  K(t,  f)  are 
related  to  the  removal  of  the  zero  and  first  degree  harmonics  from  the  gravity 
anomaly  field.  For  our  (low  altitude)  applications,  however,  these  terms  have 
negligible  effects  and  so  it  is  sufficient  to  retain  only  the  first  term  of  the 
integral  kernel,  giving: 


A  g  (  r ,  $  ,  l  ) 


tli 

4rr 


ii 


ag(R,  *’ 

D3 


2tLl 


do- 


(2.1.13) 
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3. 2  Rat io na  1  e  Beh in d  art  I nd ircct  Me thod  o f  Upward  Conti n u ation  of  Sur f ace 

Anonta  lies 

We  are  given  free-air  anomalies  Ags  on  the  earth’s  surface, and  we  want  to 
upward  continue  them.  The  difficulty  with  the  use  of  Poisson  integral  to  solve 
this  problem  is  that  the  Poisson  geometry  requires  that  the  data  to  be  upward 
continued  refer  to  points  on  the  surface  of  a  sphere.  Conceptually,  what  could 
be  done  would  be  to  first  analytically  continue  the  surface  anomalies  4gs  to 
level  surface  anomalies  Agp.  A  simple  conceptualization  of  this  continuation  is 
by  Taylor  series,  wherein  certain  corrections  that  depend  on  the  vertical 
gradients  of  the  field  are  added  to  Ags  to  arrive  at  a g|  (see  equation  (1.1)  of 
the  Introduction): 


^4 


dag 


(Hr 


Hs)  + 


1 

2  JH2 


(Hr 


(3.2. 1 ! 


It.  is  to  be  expected  that  the  rougher  the  anomaly  field  (from  which  4gs  and 
,*  or-  samples)  the  larger  will  be  the  difference  between  igs  and  igp, 
because  high  degree  frequencies  are  affected  much  by  downward/upward 
continuation  (see,  for  example,  Rummel,  1975,  pp.  42-43)  while  low  frequencies 
are  not  as  critically  affected. 

There  is  a  theoretical  problem  in  the  case  of  downward  continuation  of  Ags 
(igp  lies  below  Ags),  namely,  that  very  high  irregularities  in  the  field  caused 
by  very  high  irregularities  in  the  topography  may  cause  the  solution  igp  to 
diverge,  or,  at  least,  be  very  unstable  in  the  sense  that  small  errors  in  4gs 
will  amplify  tremendously  to  errors  in  4 gp.  Another,  now  practical,  problem  is 
that  the  computation  of  the  vertical  gradients  of  the  field,  even  just  the  first 
gradient  and  the  more  so  the  higher  the  gradients,  places  rather  severe 
requirements  on  the  density  and  accuracy  of  the  given  ags. 

The  above  problems  reduce  for  the  case  of  a  smooth  anomaly  field.  The 
corrections  (igp-4gs)  will  be  small,  the  downward  continuation  stable,  and 
simpler  computational  procedures  for  the  corrections  may  be  devised.  It  is 
then  the  central  strategy  of  what  we  would  call  the  indirect  method  of  upward 
continuation  to  explnin  away  most  of  the  high  frequencies  present  in  Ags  as 
being  the  effect  igf-  (see  si  ction  2.5)  of  certain  shallow  topographic  masses,  in 
order  to  be  left  only  with  a  relatively  low  frequency  residual  field  (Ag^-ig*) 
which  can  then  be  upward  continued  less  problematically  by  the  Poisson 
integral.  The  missing  upward  continuation  of  Ag*-  to  space  points  is  then 
carried  out  essentially  by  an  equivalent  source  technique,  which  says  that  the 
field  for  which  Ag*  are  boundary  values  has  the  topographic  masses  as 
"sources"  and,  therefore,  the  said  field  can  be  generated  by  direct  integration 
of  gravitational  influences  of  the  masses  (see  section  2.5).  In  the  next  section 
we  will  detail  the  equations  that  can  be  used  to  implement  an  indirect  method 
of  upward  continuation  of  surface  anomalies  Ags. 

It.  is  to  tie  noted  that  the  Poisson  integration,  with  its  associated  problem 
>f  requiring  data  to  be  located  on  a  level  surface,  can  be  altogether  avoided 


have  found  that  under  the  approximation  (3.1.4)  the  Ag  as  operationally 
computed  from  (3.1.1)  is  a  free-air  anomaly  on  the  earth’s  surface.  Using  a 
subscript  "s"  to  denote  a  surface  free-air  anomaly  we  finally  have  the 
interpretation: 


Ag  3  *gs(h,  ♦  ,  X) 


(3. 1.6) 


The  error  of  the  interpretation  (3.1.6)  is  given  by  the  last  term  of  (3.1.2) 
arising  from  the  difference  between  tl.3  normal  height  H*  of  P  and  the  ortho 
metric  height  H  of  P: 


Ags  h,  ♦,  X)  -  Ag  = 


dy 

a  h  J  , 


(H*  -  H) 


(3. 1.7, 


An  estimate  of  this  error  can  be  obtained  using  the  standard  normal  gradient 


J  h 


0.3086  mgal/m 


i  3 .  1  .  8 ) 


md  an  approximate  formula  for  (H*  -  II)  found  in  Heiskanen  and  Moritz  (1967, 
-.<;< ■  t.ion  8- 1 3 ) : 


'1*  I*  meters)  1gRA(gals)  '  ^(km)> 
where  a  Kit  A  ls  Bouguer  anomaly  given  by 


AgBA  :  iP.  2rkpl!. 


(3.1.10) 


A  standard  value  for  2^k  p,  where  k  =  Newtonian  gravitational  constant  arid 
p=density  of  topographic  masses,  is 


2 ekp  0 .  1119  mgal/ra , 


(3.1.11) 


corresponding  to  k=66.7  x  10~*  cm3/g/sec2  and  p=2.67  g/cm3.  Using  equations 
(3.1.8)  to  (3.1.11)  into  (3.1.7)  and  using  gravity  anomaly  and  elevation  data  in 
our  test,  area  in  New  Mexico  we  found  that  the  error  t  has  a  maximal  value  of 
0.2  nigal,  occurring  over  mountainous  terrain  of  the  area. 
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p 


vertical  gradient  of  the  normal  gravity  at  P 


?p’(0,  ♦.  *)  normal  gravity  on  the  ellipsoid,  at  the  point  P’  that  has 
the  same  (♦,  X)  as  P. 


At  this  point,  assume  that  we  do  not  know  the  vertical  location  of  the  point  at 
which  the  ag  found  by  (3.1.1)  applies,  and  let  us  find  this  location.  Let  H*  be 
the  normal  height  of  P  as  defined  in  equation  (3.1);  we  can  then  perform  the 
following  manipulations  on  (3.1.1): 


ag  -  gp(h,  X) 


♦  .  X)  + 


=  8p(h,  ♦,  X)  -  [rp'<0,  X)  +  (fj)pH  .  (^)pH*  -  (jj)pll*] 

Sp(h,  ♦,  X)  -  [xp* (0,  »,  X)  .  (jj)pII«]  *  (f£]p  (H*-H)  S3. 1.2) 


The  quantity  in  brackets  is  the  normal  gravity,  upward  continued  from  the 
ellipsoidal  point  P’  to  the  normal  point  Q  of  P  defined  in  equation  (3.1): 


rq-'H*.  ♦,  X) 


(3.1.3) 


To  a  good  approximation  we  may  assume 


H  s  H* 


(3.1.1' 


(see  equation  (3.1.7)  below  for  an  estimate  of  the  actual  difference,  H-H*). 
Substituting  (3.1.3)  into  (3.1.2)  and  neglecting  the  small  third  term  because  of 
(3.1.4)  we  get: 


ag  s  gp(h,  X)  rQ(H*>  ♦,  x) 


•'3.1.5) 


The  right  side  of  (3.1.5)  is,  according  to  (3.1),  the  gravity  anomaly  at  the  point 
P,  which  i r  the  present  case  is  a  station  on  the  earth’s  surface.  Therefore  we 


I.  Upward  Continuation  of  Surface  Free-Air  Anomalies 


In  this  section  we  will  describe  u  procedure  that  can  be  used  for  the 
upward  continuation  of  free-air  anomalies  given  on  the  surface  of  the  earth. 
These  surface  free-air  anomalies  are  boundary  values  of  free-air  anomalies  in 
space.  Recall  from  Heiskanen  and  Moritz  (1967,  pp.  91,  292)  that  the  gravity 
anomaly  on  or  above  the  surface  of  the  earth  is  defined  as  follows: 


4gp(h,  X)  -  gp(h,  <t>,  X)  -  7q(H*,  4>,  X)  (3.1) 

where  (see  diagram  on  p.  13) 

Agp  gravity  anomaly  at  P 

h  height  of  P  above  the  reference  ellipsoid 

♦,  X  geodetic  latitude  and  longitude  of  P  and  0 

Q  normal  point  of  P;  point  Q  is  established  such  that  the  actual 

gravity  potential  Wp  at  P  is  equal  to  the  normal  gravity 
potential  Uq  at  Q,  and  P  and  Q  lie  on  the  same  plumb  line  of 
the  normal  gravity  field. 

H*  height  of  the  normal  point  Q  above  the  reference  ellipsoid; 

H*  is  called  the  normal  height  of  P. 

gp  gravity  at  P 

7Q  normal  gravity  at  Q. 


Helow  we  first  show  that  operationally  available  anomalies  can  be  closely 
interpreted  as  surface  free-air  anomalies,  and  then  we  describe  a  strategy  for 
upward  continuation  of  these  anomalies  combining  Poisson  integration,  spherical 
harmonics,  and  topographic  mass  effects. 


3.J  Availnbje  Anomalies  as  Surface  Values 

In  practice,  available  gravity  anomalies  had  been  computed  from: 

Ag  --  gp( h,  ♦,  X)  -  (^)pH  -  7p>  (0,  *,  X)  (3.1.1) 

where 

h,  ♦,  X  ellipsoidal  height,  geodetic  latitude,  and  geodetic 

longitude  of  a  gravity  station  P 

gpi'h,  ♦,  X;  measured  gravity  at  P 

H  orthometric  height  of  P 
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Input:  (♦,  X,  H*)  plus  parameters  of  the  normal  ellipsoid. 
Compute: 


X  =  (N  +  H*)  cos  ♦  cos  X 

(2.6.2) 

Y  =  (N  +•  H*)  cos  ♦  sin  X 

(2.6.3) 

Z  =  (N(l-e2)  +  H*)  sin  ♦ 

(2.6.4) 

with  N  =  a/(l-e2  sin2  ♦)* 

(2.6.5) 

p2  =  X2  +  Y2 

(2.6.6) 

r2  =  p2  +  Z2 

(2.6.7) 

c  =  ae 

(2.6.8) 

K2  =  ra  +■  c2 

(2.6.9) 

h2  =  K4  -  4p2c2),t 

(2.6.10) 

.  ,  -  K2-h2 

Sln  *  '  2p2 

(2.6.11) 

tan  p  = 

p  cos  a 

(2.6.12) 

q  r  X  O  -  3  cot  a  (1 -a  cot  a)] 

(2.6.13) 

,  3(1  -  a  cot  a) 

sin2  a 

(2.6.14) 

w  =  (1  -  sin2a  cos 2P)* 

(2.6.15) 

KM  sin2«  w2a2q’  sin2<*  ,  .  _ _  1, 

r“  •  c*w  *  -  2cq.„ -  -  3 

(2.6.16) 

r  w2a2q  sin  <*  sin  P  cos  P 

P  cq0w 

(2.6.17) 

_  „  w2c  cos20 

u  a  w  tan  a 

(2.6.18) 

7p  =  -r b  +  w!c-s.in  cos  p 

r  r  w  sin  a 

(2.6.19) 

71  2  \7a  +  7p)% 

(2.6.20) 

where  yj  is  the  normal  gravity  in  the  vertical  direction. 
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geop  W  =  Wp  =  const. 


Diagram  for  Section  2.6 

The  relation  between  the  geometric  height  h,  normal  height  H*  and  the  height 
anomaly  f,  for  a  given  computation  point  P  in  space. 


For  the  computation  of  normal  gravity  y we  require  the  normal  height  H*. 
Based  on  (Heiskanen  and  Moritz,  1967,  eq.  (8-5))  H*  can  be  derived  from  h  and 
height  anomaly  <•  according  to  the  formula: 


H*  =  h  -  <•  (2.6.1) 


For  our  balloon  project  we  converted  the  known  geometric  heights  into  the 
normal  heights  using  the  height  anomalies  estimated  from  the  spherical 
harmonic  expansion  of  gravity  field  up  to  degree  180  (Rapp,  1981). 

Here  we  will  give  a  summary  of  the  procedure  as  presented  in  (Rapp, 
January  1966,  pp.  14-16): 
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anomalies  at  the  point  Q.  This  means  that  we  are  viewing  igQ4,  and  4gQ®A  as 
components  of  excesses  or  deficiencies  of  the  actual  gravity  at  the  point  Q 
over  the  normal  gravity  at  the  "corresponding"  normal  point  Q’  (the 
correspondence  between  Q  and  Q’  is  that  the  actual  potential  at  Q  is  equal  to 
the  normal  potential  at  Q’,  and  Q  and  Q’  lie  on  the  same  plumb  line  of  the 
normal  field  -  see  Heiskanen  and  Moritz,  1967,  p.  83).  In  this  sense  we  may 
also  call  agQ1  and  4gQ^A  as  gravity  anomalies  generated  by  the  topographic 
masses  and  Bouguer  plate,  respectively. 


2.6  Normal  Gravity  in  Space 

The  spatial  relationship  between  the  actual  gravity  field  and  the  normal 
gravity  field  generated  by  a  given  ellipsoid  of  reference  is  given  on  the 
diagram  on  the  next  page.  The  geop  passing  through  point  P  has  the  same 
constant  potential  as  the  spherop  passing  through  the  point  Q,  where  the  geop 
and  spherop  are  the  equipotential  surfaces  of  the  actual  and  ellipsoidal  gravity 
field  respectively. 

Normal  gravity  corresponding  to  a  given  value  of  gravity  anomaly  at  a 
fixed  location  P  in  space  is  defined  to  be  the  vertical  component  of  attraction 
generated  by  the  equipotential  ellipsoid  of  revolution  (rotating  with  the  same 
angular  velocity  <j  as  the  real  earth)  at  the  respective  point  Q  located  on  the 
equipotential  surface  of  the  ellipsoidal  field  corresponding  to  point  P.  The 
spatial  correspondence  between  P  and  Q  is  uniquely  determined  by  the 
requirement  that  the  earth’s  gravity  potential  at  P  is  equal  to  the  normal 
gravity  potential  of  ellipsoid  at  the  corresponding  point  Q.  The  normal  gravity 
in  space  is  fully  determined  by  the  geometric  (size  and  shape)  and  the 
physical  (surface  potential  and  the  rotation)  properties  of  the  level  ellipsoid. 
Combined  with  the  gravity  anomaly  the  normal  gravity  can  be  used  to  compute 
the  vertical  attraction  due  to  the  actual  Earth  at  any  location.  For  the 
purpose  of  this  report  we  use  the  equations  by  (Hirvonen,  1960),  as 
implemented  by  (Rapp,  Feb.  1966)  in  his  FORTRAN  subprogram  ’SGAMMT’. 


Although  the  method  of  Hirvonen  was  fully  described  by  (Rapp,  Jan.  1966) 
and  then  fully  documented  for  the  computer  implementation  in  (Rapp,  Feb.  1966) 
we  decided  to  state  here  the  equations  used  by  the  subprogram  ’SGAMMT’.  For 
details  the  reader  is  referred  to  (Hirvonen,  1960)  and  (Rapp,  Jan.  1966). 

Suppose  we  need  the  vertical  component  of  normal  gravity  7T  cor¬ 
responding  to  the  computation  point  P  having  the  coordinates  ♦,  X,  h,  where  h 
is  the  geometric  height  of  P  above  the  reference  ellipsoid.  Then,  following 
(Heiskanen  and  Moritz,  1967,  eq.  (8-5))  the  normal  gravity  7T  should  be 
referred  to  some  point  Q  which  is  the  'normal'  counterpoint  of  P.  Point  Q  can 
be  found  by  projecting  point  P  from  the  geop  having  potential  Wp  on  to  the 
corresponding  spherop  having  normal  potential  U=Wp  (see  diagram  below).  The 
distance  between  P  and  Q  or  the  geop-spherop  separation  is  called  the  height 
anomaly  (.  The  geometric  height  of  ’normal’  point  Q  above  reference  ellipsoid 
is  defined  to  be  the  normal  height  H*  of  the  corresponding  point  P. 
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(2.5.3) 


*4  =  “p  HI 


-  rq  cos  V'pQ  _  _ 2_ 


«Q 


tQrp 


)  dVQ 


The  last  equation  gives  the  ("topographic")  gravity  anomaly  in  space  generated 
by  topographic  masses  of  assumed  density. 

For  the  generation  of  topographic  gravity  anomalies  right  on  the  surface  of 
the  earth,  a  special  treatment  is  sometimes  convenient.  Assuming  that  the 
masses  referred  to  are  those  lying  between  the  actual  topography  and  the 
geoid  we  can  perform  the  following  split  (see  ibid.,  pp.  130-132): 


=  AggBA  _  (2.5.4) 

where 

AgQ^  vertical  attraction  at  the  surface  point  Q  generated  by  topo¬ 

graphic  masses  lying  between  the  actual  topography  and  the 
geoid. 

AgQBA  vertical  attraction  at  Q  generated  by  a  Bouguer  plate  through  Q. 

tCQ  the  well-known  gravimetric  terrain  correction  at  Q,  to  account 

for  the  difference  between  the  attraction  AgQ®A  caused  by  the 
Bouguer  plate  and  the  attraction  AgQ^  caused  by  the  actual 
topographic  masses. 

In  terms  of  formulas: 

AggBA  =  27rkpHQ  (2.5.5) 

tcQ  =  %  kpRa  |  {  dor  (2.5.6) 

<T 

in  which  (Moritz,  1966,  p.  88): 


R 

mean  earth  radius 

H 

elevation  of  integration  point 

Hq 

elevation  of  computation 

point 

<j 

unit  sphere 

d<r 

element  of  solid  angle 

to 

2R  sin  f/2 

i> 

angular  distance  between 

Q  and  dv. 

The  use  of  the  symbol  "A"  in  (2.5.4)  is  intended  to  suggest  that  in  this  report 
we  are  using  the  attractions  AgQ1  and  4gQ®A  as  components  of  gravity 
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2.5  Direct  Topographic  Mass  Effect 

Topographic  heights  are  a  much  more  readily  available  and  cheap  type  of 
information  than  gravity  anomalies  themselves,  and  can  be  effectively  utilized 
in  gravity  anomaly  interpolation  and  upward  continuation  problems.  The  idea 
is  to  subtract,  from  the  originally  given  gravity  anomaly  boundary  values  on 
the  earth’s  surface,  the  gravity  anomaly  effects  caused  by  topographic  masses 
of  assumed  density.  The  residual  anomalies  are  then  smoother  and  can  be 
much  more  easily  interpolated  and  upward  continued  than  the  original 
anomalies.  The  removed  effects  of  the  topographic  masses  are  then  added  back 
at  the  interpolation  points  or  at  the  upward  continuation  points  by  direct 
integration  of  the  gravitational  influence  of  the  masses  at  those  points. 

Topographic  masses  of  assumed  constant  density  p  (e.g.,  p- 2.67  g/cm3  is  a 
standard  density  for  land  areas)  directly  generate  gravitational  attractions  at 
points  on  or  above  the  earth’s  surface.  Considering  the  topographic  masses  as 
anomalous,  we  have  the  following  "topographic"  anomalous  potential  generated 
at  the  point  P  in  space  (Heiskanen  and  Moritz,  1967,  p.  3): 


T  t 
‘P 

=  kp 

III  <4>  dVQ 

v 

(2.5.1) 

where 

k 

Newtonian  gravitational  constant 

P 

constant  density  of  the  topographic  masses 

dVQ 

element  of  volume 

V 

volume  occupied  by  the  masses 

*Q 

(rp2+rQ2  -  2rprQ  cos  ^pQ)^,  i.e.,  the  spatial 
P  and  Q 

distance  between 

V-pQ 

angular  distance  between  P  and  Q 

rP> 

rQ 

geocentric  radius  of  P,  Q. 

We  have  the  following  fundamental  relation  between  gravity  anomaly  Ag  and 
anomalous  potential  T  in  a  spherical  approximation  (ibid.,  p.  88): 

<»r 

9 

-  T 
r 

(2.5.2) 

Substituting  (2.5.1)  into  (2.5.2),  exchanging  the  order  of  integration  .and 
differentiation,  and  performing  the  differentiations,  yield: 
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where 


(2.4.3) 


a(x, 


y) 


2 

2n 


1 _ 

(x2  +  y2  +  z2)3/2 


a(x,  y)  can  be  viewed  as  the  impulse  response  function  of  the  upward 
continuation  operation  (2.4.2). 

Now,  the  2-dimensional  Fourier  transform  of  (2.4.3)  is  defined  as 
(Robinson  E.A.,  M.T.  Silvia,  1981,  p.  224): 


A(kx,  ky) 


y)  e2*i0<xx+kyy)  ^  = 


-z 


2irV  k£  +  ki 


(2.4.4) 


for  z  >  0,  i  =  V  -1 


(see  equation  7,  p.  11  and  equation  44,  p.  56,  Erdelyi  et  al.,  1954). 

This  is  the  frequency  response  or  transfer  function  associated  with 
upward  continuation  operator.  The  function  A(kx,  ky)  is  a  spatial  frequency 
function  of  two  continuous  variables  kx,  ky  representing  the  frequencies  (in 
cycles  per  unit  length)  along  x  and  y  directions. 

The  boundary  values  f(«,  p)  can  be  transformed  to  the  frequency  domain 
by  means  of  the  Fourier  integral: 


Ffkx,  ky)  -  i  i  f(x,  y)  e-i(kx«kyy)  ^  ■ 


(2.4.5) 


Using  this  transform  the  equivalence  of  the  Dirichlet  integral  (2.4.1)  in  the 
frequency  domain  turns  out  to  be  the  multiplication  of  the  transfer  function 
(2.4.4)  and  the  Fourier  transform  of  the  data  measured  at  the  boundury 
surface: 


Fz(kx,  ky)  =  A(kx,  ky)  •  F(kx,  ky)  =  e  zW  k*+ky  •  F(kx,  ky)  (2.4.6) 


Finally,  the  desired  upward-continued  function  is  obtained  by  the  inverse 
Fourier  transform  of  the  frequency  function  Fz  (Bhattacharyya,  1967): 


fz(*,  y) 


FZ(kx, 


ky)  e"27Ti(kxX  +  kyy) 


dkxdky 


(2.4.7) 
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2.4  Fourier  Transform 


In  this  section  we  briefly  introduce  the  3-dimensional  Dirichlet  problem  for 
the  half-space  and  state  its  solution.  This  forms  the  theoretical  basis  on  which 
the  Fourier  transformation  technique  can  be  implemented  in  upward  continu¬ 
ation  problems.  Fourier  technique  offers  high  speeds  in  computation,  requiring 
no  more  assumptions  than  the  planar  Poisson's  integral  method  does.  Now  we 
will  state  the  problem  theoretically  (refer  to  later  sections  for  examples  of 
practical  implementation). 

The  Dirichlet  problem  in  three-dimensions  consists  of  the  Laplace  equation 


Ilf  iff  j.  ilf 

dX2  +  3y2  dz2 


v2f  =  0 


within  some  region  V  with  boundary  surface  S,  together  with  data  prescribed 
on  S. 

We  will  use  the  usual  planar  approximation  to  the  earth’B  surface,  where 
the  boundary  surface  S  is  the  plane  z=0  and  the  volume  of  interest  V  is  the 
half-space  z>0.  The  solution  to  this  problem  is  known  in  the  literature 
(Robinson  E.A.,  M.T.  Silvia,  1981,  p.  223)  as  the  Dirichlet  integral. 


fz(x>  y)  =  f(x,  y,  z)  =  ^  j  / 


f(a,  P )  d adp 


{ (x~ <*)  *  +  (y-p) 2  +  z2]3/2 


z  >  o 


(2.4.1) 


where  f(a,  p)  -  f(a,  p,  0)  are  boundary  values  of  f  for  z=0. 

In  this  report  the  function  f  to  be  upward  continued  is  gravity  anomaly 
function.  The  gravity  anomaly  is  a  harmonic  function  in  case  of  the  planar 
Dirichlet  problem  as  stated  above.  We  use  the  fact  that  in  our  planar  case 
vertical  gravity  satisfies  the  Laplace  equation  if  the  potential  does.  The  proof 
can  be  found  in  (Robinson  E.A.,  M.T.  Silvia,  1981,  p.  213). 


The  solution  (2.4.1)  is  identical  with  the  Poisson's  integral  for  gravity 
anomalies  introduced  in  Section  2.3.  Poisson’s  integral  can  be  viewed  as  a 
limiting  case  for  the  sphere  of  radius  R  when  R-V*.  On  the  other  hand  the 
Dirichlet  solution  (2.4,1)  is  derived  for  the  planar  case  of  half-space  z>0  using 
the  Green’s  second  identity  (Robinson  E.A.,  M.T.  Silvia,  1981,  ch.  4).  Notice 
that  (2.4.1)  represents  the  2-diraensional  convolution  integral. 

Equation  (2.4.1)  can  be  written  in  the  form: 


fz(*. 


y' 


I 

.  e» 


P)  a(x-«, 


y-p)  d«  dp 


(2.4.2) 
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Equation  (2.2.3)  was  motivated  by  the  fact  that  for  a  point  on  the  ellipsoid 
(i.e.,  r  =  rg,  H0  =  0)  (2.2.3)  becomes: 


4g(rE, 


*>■  s 


E  <"'1)  E 

n=2  m=0 


(C*  cos  mX  +  S  sin  mX)  P  (sin  ♦). 
nm  nm  nm 


(2.2.4) 


The  last  equation  is  essentially  the  equation  inverted  in  Rapp,  1981,  eq.  (2) 
(except  that  Rapp  used  the  geocentric  latitude  ?  instead  of  the  geodetic 
latitude  <>)  to  compute  his  coefficients  C^m>  S^m  from  terrestrial  data 
4g(rE>  ♦«  realizing  this  we  merely  constructed  (2.2.3)  as  an  upward 

continued  version  of  (2.2.4)  under  a  spherical  approximation.  The  rationale 
behind  (2.2.3)  does  not  hold  at  lower  degrees  (say,  n  <  36)  of  the  Rapp-1981 
field,  these  degrees  being  dominated  by  satellite  derived  coefficients  Cnm,  S^m. 
Nevertheless,  we  decided  to  use  (2.2.3)  entirely,  for  n=2  to  180*;  in  any  case, 
(2.2.3)  serves  our  above  stated  purpose  of  generating  a  self-eonsistenL,  spatial 
reference  anomaly  field  from  a  set  of  spherical  harmonic  coefficients. 


2^3  Flat-Earth  Poisson  Integral 

For  our  applications  of  the  Poisson  integral  (2.1.13a)  it  is  sufficient  to  use 
a  planar  approximation,  namely  (Rapp,  1966;  Hirvonen  and  Moritz,  1963): 


ag(r,  ♦,  X) 


dx  dy 


where 

A  fixed  integration  area 

D0  (x2  +■  y2  +  H2)* 

x  R  cos  ♦’  (X’  -  X) 

y  R(0’  -  ♦) 


(2.3.1) 


Equation  (2.3.1)  represents  a  flat-earth,  space  domain  upward  continuation 
formula  for  gravity  anomalies.  This  equation  is  indicated  to  be  valid  for  a 
distance  up  to  20"  from  the  computation  point  and  up  to  an  upward  continu¬ 
ation  distance  of  250  km  (Hirvonen  and  Moritz,  1963,  p.  71). 
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n  ina  aiinjiiHiiii'i am  i>i  i 


m  /  ~r  \  \  kM 

Tn(a,  X)  =  — 

a 


E 

m=0 


(C*  cos  mX  +  S  sin  mX )  P  (sin 
nm  nm  nra 


(2.2.1) 


where 

a 


kM 

C*  ,  S 
nm  nm 


kM  -*  kM  - 
a  nm’  a  nm 

P 

nm 


u  Mally  an  equatorial  radius 
geocentric  gravitational  constant 

fully  normalized  potential  coefficients  with  even-degree 
zonal  reference  values  subtracted 

spectrum  of  T(r,  X)  on  the  sphere  of  radius  a 
fully  normalized  Legendre  functions 


The  upward  continued  gravity  anomalies  from  spherical  harmonics  are  obtained 
by  substituting  (2.2.1)  into  (2.1.4)  with  R=a,  then  summing  the  surface 
harmonics: 


Ag'r,  ?,  X)  -  - 


kM 

2 


E 


(n 


-l)(v)n+2 


E<s 

ra=0 


cos  mX  +-  S  sin  mX)  P  (sin 
nm  nm  nm 

(2.2.2' 


For  the  purposes  of  this  report  we  would  like  to  use  spherical  harmonics 
to  generate  a  rigorously  self-consistent  field  of  gravity  anomalies  on  the 
earth's  surface  and  in  space.  This  field  is  to  be  used  as  reference,  to  be 
subtracted  from  the  observed  field  as  part  of  upward  continuation  procedures. 
For  this  purpose  (2.2.2)  can  be  used,  but  during  our  applications  and  because 
we  used  the  Rapp  (1981)  field  we  decided  to  use  the  following  equation 
instead: 


Ag(r,  ♦,  X) 


kM 


Y~  (n-l) 
n~2 


_a _ .  n  + 

a+H  } 


(C*  cos  mX  +  S  sin  mX)  P  (sin  $),  (2.2.3) 

nm  nra  nm 


where,  as  in  (2.1.13b),  H0  is  the  height  of  the  computation  point  (r,  <1>,  X)  above 
the  reference  ellipsoid. 
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Note  that  (2.1.13)  needs  gravity  anomalies  to  be  given  on  the  surface  of  a 
(geocentric)  sphere.  In  practice,  we  have  gravity  anomalies  given  on  the 
surface  of  a  (geocentric)  reference  ellipsoid  (the  problem  of  the  topography  is 
taken  care  of  separately).  To  still  use  (2.1.13)  in  practice,  we  follow  the 
spherical  approximation  used  in  Heiskanen  and  Moritz  (1967,  p.  241): 


Spherical  Approximation 


**<r.  ♦.  X)  =  |  J  X')  „„ 


t  = 


R+Hfl 


cos  ii  =  sin  ♦  sin  ♦’  +•  cos  ♦  cos  cos(X*-X) 


(2.1.13a) 


(2.1.13b) 


(2.1.13c) 


where 

(r,  ♦,  X) 

C*~E ,  ♦  X’) 

R 

Ho 


geocentric  radius,  and  geodetic  latitude  and  longitude  of 
the  computation  point  in  space. 

geocentric  radius,  and  geodetic  latitude  and  longitude  of 
data  point  on  the  reference  ellipsoid. 

a  mean  earth  radius,  taken  as  R=6371  km. 

height  of  the  computation  point  (r,  ♦,  X)  above  the  refer¬ 
ence  ellipsoid;  we  will  also  call  H0  the  upward  continua¬ 
tion  distance,  i.e.  the  distance  through  which  the  data  are 
upward  continued  to  arrive  at  the  value  of  anomaly  in 
space. 

still  evaluated  by  (2.1.11),  but  now  using  t  from  (2.1.13b) 
and  cos  i>  from  (2.1.13c). 


2.2  Sp herical  Harmonics 

The  surface  spherical  harmonics  of  the  anomalous  potential  T(r,  X)  on  a 
sphere  of  radius  "a"  can  be  written  as  (Rapp,  1982): 


with  the  use  of  either  least  squares  collocation  or  Bjerhamraar-type  equivalent 
source  technique.  In  these  methods  it  would  still  be  advisable  to  reduce  the 
original  anomalies  Agg  to  the  smoothed  anomalies  (Ags-Agt)  in  order  to  maximize 
convergence  and  economy  of  solution.  The  residual  anomalies  (Ags-Agl)  are 
essentially  inverted  into  a  new  set  of  parameters  (the  solution  vector 
in  collocation;  the  fictitious  quantities  Ag*  on  the  Bjerhammar 
sphere,  in  the  Bjerhammar  method)  and  then  the  new  parameters  rigorously 
generate  the  external  gravity  field.  The  rigor  of  the  collocation  or  Bjerhammar 
approach  lies  in  that  they  can  handle  the  fact  that  the  data  (AgB-4gl)  are 
located  on  a  non-level  surface.  The  main  disadvantage  of  these  methods  is  that 
they  require  expensive  matrix  inversion. 

3.3  Direct  and  Indirect  Upward  Continuation 

A.  Direct  Method 

In  what  we  will  call  the  direct  method,  either  the  surface  free-air  anomalies 
Ags  or  the  terrain-corrected  free-air  anomalies  (Faye  anomalies)  4gs+tc  are 
input  directly  in  the  Poisson  upward  continuation  integral  to  model  the  external 
gravity  anomaly  field.  Using  the  planar  approximation  (2.3.1)  we  have  the 
following  directly  upward  continued  fields: 

(igs)H0  =  UP{4gs>  =  2^  I  f  ^  dy  O.3. 1} 

A  ° 

and 

(Ags  +  tc)n0  =  Up((Ags  +  tc)>  =  2n  J  1  ^  dy  (3.3.2) 


where  the  superscript  D  denotes  the  direct  method;  the  subscript  H0  denotes 
the  upward  continuation  distance  H0;  Up  denotes  the  Poisson  upward  con¬ 
tinuation  operation;  Ags  is  given  by  (3.1.1);  and  tc  is  formally  given  by  (2.5.6). 
Equation  (3.3.1)  is  the  usual  simple-minded  application  of  the  classical  upward 
continuation  solution.  Equation  (3.3.2)  is  the  well-known  Pellinen  type  of 
approach  (see  details  in  the  next  paragraphs)  in  which  a  first  order  reduction 
of  surface  data  Ags  to  a  level  surface  is  implemented  using  the  terrain 
correction  and  an  assumption  of  strong  correlation  between  Ags  and  elevations 
(Moritz,  1968,  pp.  1-2). 

Note  that  in  both  equations  (3.3.1)  and  (3.3.2),  the  vertical  position  of  the 
level  surface  to  which  the  input  anomalies  are  assumed  to  refer  has  no 
clear-cut  theoretical  definition.  It  is  only  implicitly  required  that  this  level 
surface  be  close  to  the  earth’s  surface,  in  order  to  minimize  the  differences 
between  the  surface  anomalies  AgB  and  the  level  surface  anomalies  Ag*.  A 
natural  choice  for  the  position  of  the  reference  level  would  be  that  of  some 
mean  elevation  in  the  area  covered  by  the  surface  anomalies.  In  the  section  on 
numerical  investigations  (section  9)  we  present  a  study  of  the  sensitivity  of 
the  upward  continuation  results  to  the  choice  of  the  position  of  reference  level. 
Given  the  reference  level,  the  upward  continuation  distance  H0  to  be  used  in 
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(3.3.1)  and  (3.3.2)  follows  as  being  equal  to  the  elevation  of  the  computation 
point  above  the  reference  level  surface. 

The  anomalies  Ags  are  anomalies  on  the  surface  of  the  earth,  and  therefore 
the  level  surface  to  which  Ags  refers  is  strictly  changing  from  point  to  point. 
These  data  with  changing  level  are  therefore  not  directly  useable  under  the 
problem  formulation  of  the  original  Poisson  integral,  namely,  that  the  data  must 
be  on  the  surface  of  a  single  level  sphere.  However,  compared  with  Ags,  the 
quantities  (Ags  +  tc)  are  closer  to  being  level  surface  anomalies.  This  use  of 
the  terrain  correction  to  approximate  a  reduction  of  surface  data  to  a  level 
surface  is  discussed  in  Moritz,  1966,  pp.  104-107.  There  it  is  shown  that  for 
long  wavelengths  (n,  small)  we  have: 


( tc)n 


f-lig 

l  3H 


(3.3.3) 


where  AH=HS-Hp  is  the  vertical  distance  between  the  surface  anomaly  Ags  and 
the  level  surface  anomaly  Ag*  and  the  subscript  n  denotes  the  n*-*1  surface 
harmonic  of  the  quantity  in  parentheses.  We  therefore  have  this  relation 
between  the  harmonics  for  relatively  small  n  (see  (3.2.1)): 


A£n  “  (Al?s  "an  4^n  r  +  tc)n  (3.3.4) 

where  Agj^  is  the  n1-*1  harmonic  of  the  level  surface  anomalies  Ag*.  Rquation 
(3.3.4)  with  the  provision  that  n  is  small,  means  that  in  the  space  domain  the 
use  of  a  (Ags  +  tcs),  where  now  we  let  tcs  denote  a  long  wavelength  form  of 
the  actual  tc,  serves  to  implement  a  first  order  long  wuvelength  reduction  of 
the  surface  data  Ags  to  level  surface  data  Ag*.  Again,  as  stated  in  the 
previous  paragraph,  the  level  of  Ag*  is  not  clearly  defined  under  this 
“te-technique"  of  data  reduction,  and  the  sensitivity  of  upward  continuation 
results  to  a  particular  choice  of  reference  level  for  Ag*  will  be  studied  in 
section  9. 

The;  use  of  tc  instead  of  a  tcs  in  (3.3.2)  forms  a  theoretical  objection  to 
this  equation,  because  according  to  the  last  paragraph  the  interpretation  of  the 
terrain  correction  as  a  data  reduction  to  a  level  surface  is  theoretically 
guaranteed  to  be  valid  (via  (3.3.4))  only  at  long  wavelengths.  As  one  of  the 
desirable  features  of  the  indirect  method  to  be  discussed  next  this  objection  is 
minimized  because  a  tcs,  i.e.,  a  long  wavelength  form  of  tc,  is  used  instead  of 
the  tc  itself. 

B.  Indirect  Method 

In  order  to  explain  the  indirect  method  let  us  first  define  some  quantities. 

The  surface  free-air  anomaly  is  given  in  (3.3.1)  as: 
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4gs  =  g  -  fh  H  "  7 


(3.3.5) 


A  reference  free-air  anomaly  can  be  generated  from 
to  degree  Nmax  using  (2.2.3)  with  Ho=0  (see  also  (2.2.4)): 


potential  coefficients 


4gs  = 


kM 


Nmax 

e 

n=2 


(n-1) 


e 

ra=0 


(C^jj  cos  mX  +  Snu,  sin  raX)  Pnm(sin  ♦)  (3.3.6) 


the  superscript  s  denoting  spherical  harmonics. 

The  gravity  anomaly  on  the  earth’s  surface,  caused  by  masses  of  density  p 
lying  between  the  actual  topography  and  the  geoid  is  given  from  (2.5.4)  as: 


Agt1  =  2^kpH  -  tc 


(3.3.7) 


The  gravity  anomaly  on  the  surface  of  a  reference  topography,  caused  by 
masses  of  density  p  lying  between  the  reference  topography  and  the  geoid  is 
given  analogously  to  (3.3.7)  as: 


Agta  -  27TkpHs  -  tcs 


(3.3.8) 


where  Hs  is  the  orthometric  height  and  tcs  the  terrain  correction  of  the 
reference  topography.  For  our  purposes  Hs  will  come  from  a  spherical 
harmonic  expansion  of  topography  to  degree  and  order  Nmax,  corresponding  to 
the  maximum  degree  and  order  of  the  reference  anomalies  Ags  of  (3.3.6): 


^max  n 

11s  =  e  e  'nm  cos  ro* 

n=0  m=0 


Bju,,  sin  mX)  Pnm  (sin  ♦) 


(3.3.9) 


The  gravity  anomaly  caused  by  positive  and  negative  masses  of  density  p 
lying  between  the  actual  topography  H  and  the  reference  topography  Hs  can  be 
written  as: 


igt  -  Agf*  -  Agi-2 


(3.3. 10) 


Note  that  Agt'  refers  to  a  point  on  the  earth’s  topography,  whereas  Agt2  refers 
to  a  point  on  the  reference  topography.  Since  Ag*2  is  a  smooth  field  it  is 
reasonable  to  assume  that  the  analytical  continuation  of  Ag*-2  to  the  position  of 
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Agt*  js  approximately  equal  to  AgfJ  itself.  In  this  case  (3.3.10)  therefore  gives 
a  A g*-  that  refers  to  the  earth's  surface. 

Substituting  (3.3.7)  and  (3.3.8)  into  (3.3.10): 


Agt  =  -  Hs)  -  (Lc  -  tcs) 


(3.3.  IP, 


Now  a  residual  anomaly  Agr  can  be  obtained  by  removing  from  the  surface 
anomaly  Ags  of  (3.3.5)  the  spherical  harmonic  anomaly  Ags  of  (3.3.6)  and  the 
topographic  anomaly  Agf  of  (3.3.11): 


Agr  =  ~  Ags  “  4gt-  or 


(3.3. 12) 


Agr  =  Ags  -  4gs  -  2rrkp(H  -  Hs )  +  ( tc  -  tcs)  (3.3.13) 


If  the  reference  quantities  (those  with  superscript  s)  did  not  appear  in  (3.3.13) 
then  this  equation  would  give  the  expression  for  a  refined  Bouguer  anomaly 
(Heiskanen  and  Moritz,  1967,  p.  132);  but  because  of  the  presence  of  the 
reference  quantities  we  will  call  Agr  the  residual  refined  Bouguer  anomaly. 
Equation  (3.3.13)  states  that  the  original  anomalies  Ags  are  de-trended  (1.)  in 
the  long  wavelength,  by  subtracting  free-air  anomalies  Ags  generated  from 
spherical  harmonics  and  (2.)  in  the  short  wavelength,  by  doing  a  "Bouguer 
reduction”  not  with  respect  to  the  geoid  but  with  respect  to  the  higher  order 
but  still  smooth  surface  Hs  from  spherical  harmonics.  Since  Agt*  of  (3.3.12)  is  a 
smooth  quantity  (as  will  be  shown  numerically  later)  the  point  on  the  earth’s 
topography  at  which  Agr  applies  can  be  moved  vertically  to  the  point  on  the 
reference  topography  Hs,  so  that  for  subsequent  processing  Agr  is  assumed  to 
lie  on  the  reference  topography. 

Considering  the  above  definitions,  what  we  will  call  the  indirect  upward 
continuation  method  then  takes  place  as  follows: 

(1)  The  residual  refined  Bouguer  anomalies  agr  of  (3.3.13)  are  terrain 
corrected  by  tcs,  then  upward  continued  using  the  Poisson  integral  given 
by  (2.3.1): 


(Ag'"  t  tc 


.s  i  D 

H0 


Up f Agr 


tcs ! 


Hfl  f  f  (Agr  ♦  _tcs, 

2"  J  J  n3 


d.v  dy 


'3.3.  M  ) 


where  the  superscript  D  denotes  the  fact  that  (dgr  +  tcs)  is  input  directly 
into  Poisson  integral;  H„  is  the  upward  continuation  distance  and  l’p 
denotes  the  Poisson  integral  operator.  The  use  of  (igr  +  tcs)  instead  of 
simply  Agr  is  in  accordance  with  earlier  discussion  on  the  use  of  long 
wavelength  terrain  correction  (terrain  correction  indeed  has  some  long 
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wavelength  power-see  p.  74  for  a  feeling  for  magnitudes)  to  approximate 

the  reduction  of  surface  data  to  a  level  surface  (Agr  is  data  lying  on  the 

reference  topography  Hs).  Note  from  equation  (3.3.13)  that  the  use  of 
(Agr  +  tcs)  amounts  to  the  transposition  of  tcs  to  the  left  side  of  the 
equation;  this  significantly  reduces  the  computational  effort  needed  for  the 
evaluation  of  the  right  hand  side  and  is  therefore  a  definite  practical 

advantage. 

(2)  The  spherical  harmonic  anomalies  Ags  of  (3.3.6)  are  implicitly  upward 

continued  using  (2.2.3): 


=  Us{Ags} 
Ho 


""SS  (  a  V 

E  (n_1)  tt) 


(cXm  cos  mX  +  Snin  sin  mX)  Pnm(sin  ♦) 


(3.3. 15) 


where  Us  denotes  the  upward  continuation  of  the  spherical  harmonic  series. 

(3)  The  topographic  anomalies  Agt  of  (3.3.11)  are  implicitly  upward  continued 
by  integration  of  the  gravitational  attractions  caused  by  the  masses 
generating  Agl,  namely,  the  masses  lying  between  the  actual  topography 
and  the  reference  topography  (see  equation  (2.5.3)).  In  practice  the 
integration  can  be  implemented  using  prisms  as  integration  elements.  This 
prism  integration  is  implemented  in  an  operational  program  by  Forsberg 
(see  section  7).  In  symbolic  form, 


Ag*.  -  UflfAg1} 
Ho 


(3.3. IB! 


where  Up  denotes  the  upward  continuation  by  prism  integration  of  masses. 

(4)  Adding  tcs  on  both  sides  of  (3.3.12),  applying  the  Poisson  integral  operator 
Up,  and  transposing  terms,  we  arrive  at: 


Up  figs  *  tcs }  --  Up(4gr  4  tcs}  i  Up(ags}  t-  Upiag1 


3.3. 17) 


The  indirect  upward  continuation  method  consists  of  replacing  the  last  two 
terms  in  (3.3.17)  as  follows: 

A  gs  4  tc*)*  U  y  f  Ags  4  tcs}  -  Up|4gr  4  tcs}  4  Us(Ags}  4  UpUgM  (3.3.  IB': 


that  is,  the  operator  Up  operating  on  anomaly  data  4gs  and  4gt  have  been 
replaced  by  operators  Ua  operating  on  spherical  harmonic  coefficients  and 
Ur  operating  on  topographic  masses.  The  superscript  I  in  (3.3.18)  denotes 
the  indirect  upward  continuation  method,  and  the  U[  symbolizes  the 
(implicit)  indirect  upward  continuation  operator.  The  last  two  terms  of 
(3.3.18)  are  given  in  (3.3.15)  and  (3.3.16). 


Starting  from  Section  7  and  onwards  we  give  the  relevant  numerical 
studies  on  the  direct  and  indirect  upward  continuation  methods,  as  well  as 
studies  on  the  Fourier  transform  method  of  upward  continuation.  Meanwhile,  in 
the  next  three  section  (4,  5,  6)  we  present  a  general  global  study  of  truncation 
theory  for  the  anomaly  field,  spectral  characteristics  of  the  anomaly  signal,  and 
error  analysis  for  the  (Poisson)  upward  continuation  operation. 


4.  Anomaly  Field  Truncation  Theory  at  Altitude 

We  can  collect  equations  from  section  (2.1)  as  follows: 


Agr(r,  x)  =  ^  J  J  K(t,  *)  Ag(R,  X')  d* 


(4.1) 


K(t,  f)  =  -  t3  -  3t3  cos 


(4.2) 


K(t,  f)  =  E  (2n+l)  tn+3  Pn(cos  i>) 
n=2 


(4.3) 


A«r(r,  X)  =  tn+a  Agn(R,  ?,  X);  t  =  (f) 

n=2 


(4.4) 


A  "truncated"  gravity  anomaly  field  can  be  generated  as  follows: 


*Ir  =  4“  J  J  K(t,  i>)  Ag  d<r 


(4.5) 


K(t,  f) 


0  *  1>  *  1 0 


K(t ,  i>)  ,  i>0  *  i>  *  n 


(4.6) 


T^o  •••  truncation  cap  radius. 


This  truncated  field  is  generated  by  data  function  values  Ag  outside  a  cap  of 
radius  i>0  centered  at  the  computation  point.  In  this  sense,  this  field  is  really 
a  "remote  zone"  field  being  generated  by  remote  zone  data.  The  truncation 
kernel  can  be  expanded  into  a  series  of  Legendre  polynomials  as  follows: 


K(t,  f)  =  ^  {2n+L)  Pn(cos  f) 

n=0 


(4.7) 
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and  the  frequency  domain  generation  of  the  truncated  field  is  then  given  by 
(6ee  analogous  manipulations  in  Jekeli,  1980,  for  example): 


**Kr,  ♦.  *)  =  |  YZ  ^  Agn(R’  ?*  X)  (4'8) 

n=2 

The  truncation  coefficients  qn  can  be  obtained  as: 

qn  =  |  K(t,  V-)  Pn(cos  i>)  sin  V-  d 1>  (4.9) 

O 


=  |  K(t,  i>)  Pn(cos  i>)  sin  f  <H>  (4.10) 

■to 


Putting  y  =  cos  1>,  yc  =  cos  ^o.  we  have: 


9n( 


Yo 


K(t,  y)  Pn(y)  dy 


(4.11) 


Substituting  (4.2)  into  (4.11)  and  using  recursive  integral  evaluations  found  in 
Sheppard  (1979,  p.  E-l)  we  arrive  at  the  following  recursive  computations  for 
the  truncation  coefficients: 


qn(t,  y0)  =  tJ(l-tJ)  LnU,  y0)  -  ta  In(y0)  -  3t3  Hn(y0) 


(4.12) 


Ln(t.  yo)  =  J  pn(y)  dy 


(4.12) 


-1 


f  y0 

in(fc»  yo)  =  I  pn(y)  ^ 
-l 


(4.14) 


f  Yo 

Hn( t ,  y0)  =  J  y  pn(y)  dy 


(4.15) 


The  above  recursive  computations  were  checked  for  correctness  and 
stability  against  other  published  results  (with  excellent  agreement)  as  part  of 
the  studies  of  Cruz  (1985)  on  truncation  coefficients  for  various  gravimetric 
quantities  in  space.  The  recursive  formulas  (4.18)  for  Ln(t,  yQ)  which  is 
probably  the  main  source  of  instability  of  the  recursions  can  be  derived  either 
analogously  to  the  way  Shepperd  (1979,  pp.  B-l  to  B-3)  derived  his  Kn(t,  y0) 
functions,  or  as  a  special  case  of  a  general  formula  given  by  Jekcli  (1982, 
equations  (18)  to  (22)). 

The  truncation  method  expressed  by  equations  (4.5)  and  (4.6),  when;  the 
original  kernel  K ( t,  f)  is  set  to  zero  for  0  *  y  *  f0  is  in  accordance  with  what 
is  referred  to  as  an  unmodified  Molodensky  truncation  method.  Other  methods 
of  truncation  can  be  defined  by  specifying  the  truncation  kernel  as  in  (4.6), 
and  deriving,  using  (4.9),  the  truncation  coefficients  that  will  enter  (4.8). 
Well-known  alternative  methods  for  generating  a  truncated  field  include  the 
Meissl  truncation  method  and  an  improved  Molodensky  truncation  method,  all 
these  being  discussed  in  Jekeli  (1982). 
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The  interest  in  truncated  fields  usually  lies  in  trying  to  generate  the 
various  truncated  fields  (unmodified,  Meissl,  or  improved  Molodensky  truncated 
field)  from  a  finite  set  of  available  spherical  harmonic  coefficients  of  the 
earth’s  gravitational  potential.  This  is  in  contrast  to  generating  the  field  from 
a  set  of  Ag  values  on  a  geocentric  sphere.  The  spherical  harmonic  generation 
of  the  truncated  field  can  be  accomplished  by  using  in  (4.8)  the  surface  har¬ 
monics  Agn  generated  from  potential  coefficients  to  degree  Nmax  (see  Rapp, 
1982,  p.  4).  Specifically,  we  can  collect  the  necessary  equations  for  the 
spherical  harmonic  generation  of  a  truncated  anomaly  field  at  altitude,  as 
follows: 

___  .  Nmax 

*gf(r,  ¥,  X)  =  |  qn(t>  Vo)  *  4gn<R*  X>  (4.19a) 

n=2 


4gf(r*  *>  x) 


kM 

2aa 


^max 

y  4n(t,  y0)  * 

n=2 


(n-1) 


*  y  (c^ 

m=0 


cos  mX  +  Sjuj,  sin  mX) 


Pnmisin  ?) 


(4.19b) 


t 


R 

R+H0 


(4.19c) 


y0  =  cos  i>0 


(4. 19d) 


where 

r,  ?,  X 

qn(t.  yo) 
^gn(^»  ^ >  x) 

R 


geocentric  radius,  geocentric  latitude  and 
"geocentric"  longitude  of  the  computation  point; 

truncation  coefficients; 

surface  harmonics  of  the  spatial  anomaly  field,  on  a 
sphere  of  radius  R; 

a  mean  earth  radius,  taken  as  R=6371  km; 
an  average  value  of  gravity; 


run 


a 


fully  normalized  potential  coefficients  with  even 
degree  zonal  reference  values  subtracted; 
radius  of  the  sphere  to  which  and  refer, 
usually  an  equatorial  radius; 
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w 


~  1 


*I"P.  I"  »’ 


H0  height  of  computation  point  above  the  sphere  of  radius 

R; 

Vj  angular  radius  of  the  truncation  cap 

The  truncated  spherical  harmonic  field  TTg,- >  given  in  (1.19b)  contains  two 
types  of  errors:  the  commission  error,  which  is  the  error  due  to  errors  in  the 
potential  coefficients  being  used,  and  the  omission  error,  which  is  the  error- 
due  to  the  use  of  only  a  finite  set  of  potential  coefficients  and  the  consequent 
omission  of  the  higher  frequencies  which  may  exist  in  the  actual  truncated 
field  Xgr  being  represented.  The  combined  global  mean  square  value  of 
commission  and  omission  errors  can  be  expressed  as  (Jekeli,  1980,  p.  16;  see 
also  (4.19a): 

_  .  fynax  .  <° 

-  <\  qnJ  4Cn  +  ~  ^  lln  (4.201 

n  "  2  n  ~  ^max  1 


where 


9  n 


<5C 


n 


C 


n 


truncation  coefficients  as  found  from  equations  ( 4 . 1 11 )  (1 . 18); 


anomaly  error  degree  variances  referred  to  a  sphere  of  radius 
R  (R=6371  km,  the  mean  earth  radius).  ^Cn_is  caused 
by  potential  coefficient  errors  »C*ra  and  tSnm.  Considering 
(4.19b)  we  have: 


<5Cr 


IV 


,  2n+4 


l  i 

E 

mr0 


[(*C* 


nm  > 


'•Sum  ■  ‘ . 


1 4 .21a) 


modeled  anomaly  degree  variances  referred  to  a  sphere  of  radius 
R.  According  to  the  Tschern ing/Rapp  (1974)  model  (with 
R-637 1  km) : 


C 


n 


425.28  (n-1) 
( ii- 2)  (n+24) 


(0. 999617)n+2 


(4.21b) 


Optimization  of  -truncation  method  usually  means  modifying  the  definition  (4.6) 
for  the  kernel  K ( t ,  f)  (and  therefore  modifying  the  truncation  coefficients  qn) 
such  that  the  mean  square  error  m?  in  (4.20)  is  reduced  compared  with  t  he 
unmodified  case,  for  a  given  number  and  accuracy  (Nmax  and  <5Cn)  of  potential 
coefficients  and  given  model  (Cn)  for  anomaly  degree  variances.  Jekeli  (1982) 
has  shown  that  for  the  upward  continued  gravity  anomaly  field,  optimization  of 
truncation  method  by  using  the  Meissl  or  improved  Molodensky  technique  does 
not  offer  a  significant  improvement  over  the  unmodified  Molodensky  method 
defined  by  (4.5)  and  (4.6).  In  this  report  we  limit  ourselves  to  the  use  of  the 
unmodified  truncation  method. 


> 
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As  an  application  of  the  important  equations  (4.19)  and  (4.20)  of  this 
section,  we  considered  the  main  interest  in  this  report,  namely,  the  upward 
continuation  of  anomalies  in  the  New  Mexico  area  to  an  approximate  altitude  of 
28.5  km.  We  first  applied  (4.19)  by  modifying  program  F388  at  OSU  (see 
Section  7)  to  introduce  the  truncation  coefficients  qn(t,  v0).  We  used: 

•  the  Rapp-180  field  (Rapp,  1981) 

•  H0  -  28.5  krn  (upward  continuation  distance) 

•  y  -  979770  mgals  (mean  value  of  gravity) 

•  Vo  -  3*  (truncation  cap  radius) 

•  a  -  6378137  m  (equatorial  radius) 

•  R  -  6371000  m  (assumed  ground  level  radius) 


The  coefficients  qn  were  computed  by  a  subroutine  modified  from  a  subroutine 
published  by  Shepperd  (1979).  The  resulting  truncated  anomaly  field  for  the 
area  of  interest  in  our  balloon-borne  gravity  experiment  (Section  8)  is  given  as 
Figure _ 1. 

If  we  were  using  the  Rapp-180  field  to  account  for  the  remote  zone  outside 
a  3* -cap  centered  at  the  computation  point,  the  values  on  the  map  would  be 
added  to  the  result  of  integration  of  data  inside  the  3*-cap.  Note  that  in  this 
case  the  data  cap  of  integration  moves  from  computation  point  to  computation 
point.  For  our  final  operational  procedures,  however,  we  simply  used  a  fixed 
data  cap  to  compute  our  anomalies  in  space  and  neglected  remote  zone  effects: 
this  neglection  is  justified  by  selecting  a  sufficiently  large  size  for  the  fixed 
integration  cap.  In  the  case  just  considered  an  integration  cap  of  radius  3*  is 
appropriate  since  the  remote  zone  effects  are  small  as  shown  in  Figure  1. 

Let  us  now  turn  to  an  application  of  equation  (4.20).  Application  of  (4.20) 
yields  u  global  analysis  of  the  effect  of  remote  zone  anomalies  on  the  upward 
continued  anomaly  and  gives  an  indication  of  the  ability  of  spherical  harmonic 
model  to  account  for  this  effect.  As  a  first  application  we  used  only  the 
second  term  of  (4.20)  and  started  the  summation  from  n=2.  This  is  equivalent 
to  not.  using  any  spherical  harmonic  model  to  account  for  the  remote  zone 
effects,  the  total  error  being  one  of  pure  omission.  Using  the  Tscherning/Rupp 
(1974)  Cn-model  (4.21b)  the  results  of  summing  the  second  term  of  (4.20)  for 
altitude  30  km  and  various  truncation  angles  are  shown  in  Figure  2  as  the  top 
curve.  This  curve  shows  that  a  cap  of  radius  3*  is  needed  to  reduce  remote 
zone  effects  to  submilligal  level;  the  3*  cap  radius  (about  300  km  at  ground 
level)  is  about  ten  times  the  upward  continuation  distance  of  30  km  and  is 
therefore  in  accordance  with  the  rule  that  data  out  to  a  distance  ten  times  the 
upward  continuation  distance  should  be  used  (Hirvonen  and  Moritz,  1963,  p. 
68). 


LRT ! TUDE 


LONGITUDE 


Figure  1.  Rapp-180  Truncated  Gravity  Anomaly  Field  at  28.5  km  Altitude,  New 
Mexico  Balloon  Gravity  Test  Site.  C.I.  =  0.01  mgal;  Truncation  Cap  Radius  =  3". 


Figure  2.  Krror  Incurred  When  Not  Using  and  When  Using  a  Spherical 
Harmonic  Model  to  Represent,  at  Altitude  30  km,  the  Remote  Zone  Field 
Generated  by  Data  Outside  a  Cap  of  Given  Radius.  Cn-model  from 
Tscherning/Rapp  (1974);  <5Cn  from  Rapp-180  field  with  m^g  =  10  mgals  (Rapp, 
1981,  eq.  30). 


With  the  numerical  values  of  C0  =  683  (mgal2),  d  =  40  km  and  A  implied  by  5’x5’ 
blocks  we  compute 


mji  =  11.6/H2  mgal  (H  in  km). 
So  for  flight  elevation  H  =  30  km 


in„m  0.013  mgal 
30km 


which  is  a  much  smaller  contribution  than  any  of  the  effects  given  in  Table  1 
or  2.  A  spherical-earth  analysis  of  the  representation  error  for  gravity 
anomaly  upward  continuation  may  also  be  done  based  on  Sunkel  (1981,  p.  17); 
however,  the  effect  of  representation  error  is  so  small  (using  5’x5’  mean 
values)  that  we  do  not  repeat  this  type  of  analysis  here. 
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C(s)  =  M(Agag’)  is  the  covariance  function  of  actual  gravity  anomalies, 

that  is  a  statistical  characterization  of  a g  field  itself  (and  not 
the  errors  in  ig) 

'ter  (Moritz,  1962}  assume  the  following  model  covariance  function  for  gravity 
lomalies: 


C(s)  = 


lere  C0  is  the  assumed 

'  we  define  a  constant 
>.3 .5)  reduces  to: 


C( s )  C0  -  a3 


gravity  anomaly  variance. 

a  such  that  a2  -  — 9-  then  in 
d 


(6.3.5) 


first  approximation 


(6.3.6) 


isume  also  all  blocks  are  squares  having  the  area  a?  =  a2  =  A. 

hen  (Moritz,  1962)  gives  the  following  closed  expression  for  the  mean  square 
^presentation  error  (6.4.3) 


(6.3.7) 


here 


Yo 


,  9 


2 
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U3  + 
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2 
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quation  (6.3.7)  gives  the  effect  due  to  a  single  compartment.  Next,  neglecting 
gain  correlations  between  compartments  we  can  find  variance  of  the  total 
f fee t  just  by  summing  variances  of  euch  contributors: 


mK 


r  mu 

his  is 
r  ror . 


g2  A1 
384-1  U* 


(6.3.8) 


J . 0288  <«  Y  u 1  b 1  f(2  where  a  and  b  form  the  sides  of  the  block. 

the  square  root  of  error  variance  of  upward  continued  representation 
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where  Dt  is  defined  by  the  relation: 


A'  til*  ",th  01  -  1(1  '  <*-*’)*  .  (yy’> 

a, 


and  Aj  is  the  area  of  the  rectangular  block  which  is  represented  by  the  mean 
value  Ag~I, 

Define  the  error  of  representation  £  to  be  the  function  such  that  over  each 
rectangular  block  i  the  following  relation  holds: 


*  j  =  4g| .  -  4Si 

'  i  '  l 


(6.3.2) 


which  is  the  difference  between  Ag  and  the  constant  Sg[  over  the  rectangular 
block  i.  Then,  from  (6.1.2)  the  propagated  error  of  representation  due  to 
single  block  i  is: 
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terrestrial  data  should  be  avoided  as  much  as  possible,  because  these  cause 
errors  to  show  up  significantly  in  the  upward  continued  anomaly  field.  Note 
from  Table  2  that  the  larger  the  error  correlation  length  the  slower  the  rate;  at 
which  the  error  loses  its  total  energy  (variance)  with  altitude,  an  expected 

result.  I 


Table  2 

Upward  Continued  Error  Variance  and  Correlation  Length, 
for  Attenuated  White  Noise  Error  Model. 
vr  =  100  mgal2 


H=30  km 

H=10  km 

H-f 

km 

£r 

«r 

'/"v 7 

f  r 

(r 

/“vv 

(km) 

(km) 

(mgal) 

(km) 

(mgal) 

( km) 

( mgal ) 
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0.42 

17 

1.18 

10 

2.  10 

5 

50 

1.00 

20 

2.50 

12 

4.00 

10 

55 

1.82 

25 

4.00 

18 

5.71 

30 

75 

4.00 

45 

6.67 

38 

8.00 

6.3  Flat-Earth  Representation  Error  Propagation 

The  effect  of  representation  error  can  also  be  evaluated.  This  error  is 
committed  when  the  continuous  function  ag  is  replaced  by  the  step  function 
composed  of  the  mean  values  representing  it  over  rectangular  blocks  used  in 
computation.  See  (Sunkel,  1981).  Following  (Moritz,  1962)  where  the  same 
error  was  considered  (but  under  the  name  of  integration  error)  we  pose  the 
problem  in  the  following  way: 

Instead  of  the  exact  Poisson’s  integral  (6.1.1)  in  the  actual  computation  we  use 
the  summation: 


i>'.ll 


li 


A 


(6.3.  1  ) 
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from  (6.2.4).  In  Figure  5  we  show  the  covariance  function  (6.2.4)  for  r=R  (i.e., 
on  the  surface  of  the  earth)  and  varying  depths  (R-Rb)  to  the  white  noise 
process.  Figure  5  shows  clearly  the  relation  (6.2.5);  the  graphs  were  scaled  to 
yield  the  variance  vr=100  mgal2. 

The  functions  (6.2.4)  such  as  shown  in  Figure  5  can  be  used  as  error 
models,  on  the  earth  sphere,  and  it  is  of  interest  to  see  how  these  functions 
propagate  through  the  upward  continuation  operator  for  gravity  anomalies.  We 
have  the  following  procedure  for  applications: 

1.  we  have,  given,  the  error  variance  vp  and  error  correlation  length  fp  of 
the  error  process  that  we  want  to  model  on  the  earth  sphere  R. 

2.  the  depth  D  (not  to  be  confused  with  the  same  symbol  in  (6.2.2))  of  the 
white  noise  process  generating  fp  is  then  found  from  (6.2.5): 


D  =  (R-Rb)  =  |  fR 


(6.2.7) 


3.  the  constant  of  proportionality  in  (6.2.6)  is  found  from 


const.  =  vp(R-Rb)2  =  vp  D2  (6.2.8) 


4.  the  upward  continued  error  variance  at  altitude  H  above  the  earth  sphere 
R  is  found  by  applying  (6.2.6)  with  r=R+H  and  Rb=R-D: 


(6.2.9) 


5.  the  upward  continued  correlation  length  fr  is  found  from  (6.2.5): 


const . 

( r-Rb) 2 


=  VR 


D 

H+D 


(r-Rb)  -  2  (H+D) 


(6.2.10) 


We  applied  (6.2.7)  -  (6.2.10),  starting  with  vp=100  mgal2  and  various 

correlation  lengths  fp.  The  upward  continued  values  vr  and  (r  for  various 
upward  continuation  distances  are  shown  in  Table  2.  The  values  (vr)*  in  the 
table  are  directly  comparable  with  the  values  mp  in  Table  1,  and  we  see  a 
reasonable!  agreement;  note  that  (vr)  '  has  no  specific  problems  associated  with 
large  correlation  lengths  (p  whereas  mp  has  problems  with  this  as  mentioned  in 
.Section  6.1.  The  conclusion  from  Table  2  is  that  error  correlations  in  the 
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D(P,  Q)  =  v02  ^  (2n  +  1)  Pn ( cos  i>)  (6.2.2) 

t\-2 

where  we  start  the  summation  from  n=2  in  anticipation  of  inputting  D(P,  Q)  to 
an  upward  continuation  operation  that  starts  with  n=2.  D(P,  Q)  represents  a 
dirac  delta  function  (Rummel,  1982,  p.  30): 


D(P,  Q)  = 


for  P  =  Q 


0  for  P  f  Q 


P,  Q  t  a 


(6.2.3) 


that  is,  the  total  variance  is  infinite  and  the  correlation  length  is  undefined. 
To  obtain  a  finite  variance  and  correlation  length  we  upward  continue  the  noise 
process,  resulting  in  the  attenuated  covariance  function  (note  that  we  are 
using  the  gravity  anomaly  upward  continuation  operation): 


Dr(s,  VO  -  <r0 2  y  sn+2  (2n+l)  Pn(cos  V) 
n=2 


(6.2.4) 


where  s  =  (^j  ,  and 


we  visualize  that  the  white  noise  process  is  located  on  a  sphere  of  radius 
internal  to  the  earth  sphere  of  radius  R,  and  the  attenuated  white  noise 
process  described  by  (6.2.4)  is  located  on  a  sphere  of  radius  r  with  r>Rb- 

Two  important  quantities  to  characterize  the  covariance  function  Dr  are  its 
variance  vr  and  correlation  length  £r.  A  numerical  study  of  (6.2.4)  reveals  the 
following  good  approximations: 


(r-Rb)  r  3  *r 


(6.2.5) 


const . 
Vr  '  (r-Rb)2 


(6.2.6) 


that  is,  the  correlation  length  £r  is  1.5  times  the  upward  continuation  distance 
(r-Rb),  and  the  variance  of  the  attenuated  white  noise  process  attenuates  with 
the  square  of  the  upward  continuation  distance.  The  same  relations  were 
found  by  Sunkel  (1981,  p.  12,  14)  for  a  slightly  different  covariance  function 
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Now  we  can  evaluate  (6.1.6)  using  different  trial  correlation  lengths  {  at 
different  elevations  H.  The  results  are  shown  in  Table  1  using  the  parameter 
<7o  =  (10  mgal)2.  The  closed  formula  (6.1.6)  was  derived  under  the  assumption 
that  the  ratio  ( /H  remains  small  (see  Moritz,  1962,  p.  7).  For  that  reason  the 
numbers  in  the  lower  right-hand  portion  of  Table  1  have  been  crossed-out  as 
they  are  considered  not  to  have  any  physical  meaning. 

To  overcome  this  assumption  see  below  for  a  spherical  earth  analysis. 


Table  1 

Square  root  of  upward  continued  error  variance  at  different  elevations 

mn  in  [mgal]:  <r0  -  100  mgal  2 
(crossed-out  values  nave  no  physical  meaning) 


Assumed  correlation 
distance  in  error 
function  at  zero 
level  (£  in  [km]) 

H  =  30  km 

H  =  10  km 

H  =  5  km 

2 

0.28 

0.85 

1.70 

5 

0.71 

2.12 

4.25 

10 

1.42 

4.25 

20 

4.25 

.  _  . ... 

6.2  Spherical-Earth  Data  Error  Propagation 

Let  us  now  study  the  propagation  of  data  error  through  the  upward 
continuation  operation,  using  as  error  model  an  attenuated  white  noise  process 
(see  Heller  and  Jordan,  1979,  for  interesting  geodetic  applications  of  such 
process).  The  degree  variances  of  white  noise  can  be  written  as: 


dn  -  v02  (2n+l) 


(6.2.1) 


where  the  unit  variance  cr0 1  is  a  constant  and  is  equal  to  the  variance  of  a 
single  harmonic  of  degree  n  and  order  m  of  the  white  noise  process.  White 
noise  is  useful  for  approximating  an  uncorrelated  noise  process. 

The  covariance  function  corresponding  to  (6.2.1)  is: 
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1981).  Equation  (6.1.3)  shows  that  the  weakly  correlated  (white)  noise  in  the 
data  will  attenuate  very  rapidly  with  the  elevation  H,  whereas  the  long 
wavelength  widely  correlated  components  of  the  noise  will  propagate  almost 
unattenuated  into  upward  continued  gravity  field. 

This  discussion  shows  that  the  correlation  length  of  the  errors  present  in 
the  gravity  material  will  play  the  major  role  in  the  upward  continuation  error" 
analysis. 

As  a  statistically  appropriate  measure  of  the  upward  continued  error 
function  (6.1.2)  we  choose  (after  Moritz,  H.,  1962)  the  mean  square  error 

defined  as: 


"<■*«>  J|  u  •<*.  y,  *•.  y)  ^  ^ 


x  y  x  y 


(6.1.4) 


where  <r(x,  y,  x’,  y’)  =  M(e,  e’)  is  the  error  covariance  function,  which  is  a 
statistical  description  of  the  errors  s(x,  y)  in  Ag.  M  is  the  suitable  averaging 
operator. 

In  (Moritz,  1962)  it  is  shown  that  introducing  some  specific  model  of  the 
error  covariance  function  we  can  produce  explicit  expressions  for  the  upward 
continued  mean  square  error  (6.1.4). 

Suppose  we  model  the  error  covariance  function  present  in  Ag  to  be 
(Moritz,  1962,  p.  3): 


a(x,  y,  x\  y’)  =  <r0  e~c2s2  =  ff(s) 


(6.1.5) 


which  is  a  function  of  distance  only  (stationarity  and  isotropy).  In  eq. 

(6.1.5): 

S  =  -J  (x-x1)2  +  (y~y’)2  is  the  Cartesian  distance  between  two  locations 

on  the  plane. 

is  the  error  variance  (square  of  the  standard 
error  present  in  the  datum  surface  anomalies  Ag) 

The  constant  c  =  '/<n2/ {  is  inversely  proportional  to  the  correlation 
length  i  of  the  error  function  present  in  Ag  ((  has  to  be  in  the  same  units  as 
s  and  H) . 

Using  this  model  (Moritz,  1962)  shows  that  the  mean  square  error  (6.1.4) 
present  in  the  upward  continued  field  takes  the  simple  form: 


m2n 


L.  CTo 

8H2  c3 


>6.1.6) 
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Error  Propagation  for  Upward  Continuation  Operator 


6.1  Fla t-Ea rth  Data  Er ror  Propagation 

Based  on  (Moritz,  H.,  1962)  we  can  formulate  the  problem  in  the  following 


The  upward  continuation  integral  (in  planar  form)  for  gravity  anomalies  is: 


AgH(x.  y)  =  ^  J  J  4g(x\  y’) 


(6.1.1) 


where  x’y’  are  variables  of  integration, 


D2  =  H2  +  (x-x*)2  +  (y-y’): 


H  is  the  upward  continued  distance. 


Then  formally  any  error  e(x,  y)  in  terrestrial  4g  will  propagate  as: 


/  \  H  f  f  ,  ,  , .  dx’c 

:»(x-  y)  =  S  J  J  e<x  ’  y  }  “j ^ 


(6.1.2) 


Equation  (6.1.2)  has  the  exact  form  of  the  original  Poisson’s  integral  with 
gravity  anomalies  Ag  replaced  by  the  error  function  c.  As  the  computational 
point  (x,  y)  sweeps  the  particular  level  plane  at  the  elevation  H  above  the 
reference  datum  plane,  the  function  c}j(x,  y)  describes  the  variation  of  the 
directly  propagated  (upward  continued)  error. 

Formula  (6.1.2)  describes  the  sensitivity  of  the  upward  continuation 
operation  to  the  uncertainties  in  the  data.  It  assures  that  the  errors  attenuate 
according  to  exactly  the  same  law  (upward  continuation  law)  as  the  original 
data.  The  frequency  domain  equivalent  of  formula  (6.1.2)  is:  (see  eq.  2.4.4) 


EH(kx,  ky)  =  e  H2n  k*  *  ky  •  E(kx ,  ky) 


(6.1.3) 


where  > ,  E  and  tg,  are  Fourier  transform  pairs  and  kx,  ky  are  spatial 

frequencies  along  x  and  y  directions  (Robinson  E.A.,  M.T.  Silvia,  Chap.  2.4, 
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Figure  4.  Anomaly  Degree  Variances  at  Altitudes  30,  60,  100,  and  150  km. 
Tscherning/Rapp  (1974)  Anomaly  Degree  Variance  Model  Used. 
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Figure  5.  Attenuated  White  Noise  Covariance  Functions  (scaled),  for  White  Noise 
Depth  D=(2/3)(10,  20,  30,  40,  50,  100  km).  Note  that  the  Resulting  Correlation 
Length  £  =  3/2  D. 
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Figure  3a.  Gravity  Anomaly  Information  Beyond  Degree  N,  for  Altitudes  30,  60, 
100,  and  150  km.  Tscherning/Rapp  (1974)  Anomaly  Degree  Variance  Model 
Used. 


Figure  3b.  Gravity  Anomaly  Information  Generated  by  (Remote  Zone)  Data 
Outside  a  Cap  of  Given  Radius,  for  altitudes  30,  60,  100,  and  150  km. 

Tscherning/Rapp  (1974)  Anomaly  Degree  Variance  Model  Used. 


The  sum  (5.4)  is  important  for  studying  the  remote  zone  effects  on  upward 
continued  anomalies. 

Formulas  (5.1),  (5.3),  and  (5.4)  were  used  at  altitudes  30,  60,  100  and  150 
km,  the  results  being  given  as  Figures  3a,  3b,  and  4.  The  Tscherning/Rapp 
(1974)  CnR-model  was  used  (see  (4.21b)).  Figure  3a  gives  the  gravity  anomaly 
information  beyond  degree  N:  for  example,  using  a  180-field  to  represent  the 
external  anomaly  field  leaves  5  mgals  to  be  resolved  at  30  km,  1.7  mgals  at  60 
km,  0.4  mgal  at  100  km,  and  0.1  mgal  at  150  km.  Figure  3b  gives  the  gravity 
anomaly  information  due  to  data  outside  a  cap  of  given  radius:  for  example, 
using  a  truncation  radius  of  3*  leaves  a  0.9  mgal  remote  zone  anomaly  at 
altitude  30  km,  1.8  mgals  at  60  km,  2.8  mgals  at  100  km,  and  3.8  mgals  at  150 
km.  A  combined  interpretation  of  Figures  3a  and  3b  says  that  the  higher  the 
point  of  upward  continuation  the  larger  the  data  cap  needed  to  maintain  a 
desired  accuracy  level  (Figure  3b),  however,  at  the  same  time  the  resolution 
needed  for  data  inside  the  cap  becomes  less  and  less  with  altitude  (Figure 
3a)  -  this  conclusion  is  well-known  and  is  true  for  the  computation  of  all 
gravimetric  quantities  in  space.  Finally,  Figure  4  gives  the  degree  variances 
of  the  anomaly  field  at  altitudes  30,  60,  100  and  150  km  to  Bhow  the  effect  of 
upward  continuation  on  component  powers  of  the  field. 
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5.  Global  Characterization  of  Anomaly  Fields 

A  most  common  statistical  characterization  of  anomaly  fields  can  be  made  by 
giving  the  degree  variances  of  the  field,  which  gives  the  break  down  of  field 
power  by  wavelength.  Given  the  degree  variances  CnR  on  a  sphere  of  radius 
R,  the  degree  variances  on  an  external  sphere  of  radius  r  is  given  by  using 
(2.1.5)  and  analogous  derivations  to  Heiskanen  and  Moritz  (1967,  p.  260-261): 


Cnr  -  sn+J  CnR  ;  s  -  [r]  (5.1) 


The  factor  sn+2  indicates  how  each  component  wavelength  power  of  the  field 
attenuates  with  the  upward  continuation  distance  (r-R).  Note  that  the  degree 
variances  Cnp  correspond  to  the  power  of  field  features  on  a  sphere  of  radius 
R,  these  features  having  a  minimum  half-wavelength  of  approximately: 


(linear  units) 


(5.2) 


A  second  field  characterization  is  obtained  by  Bumming  the  degree 
variances  of  the  field  above  a  specified  degree  N.  This  sum  gives  a  measure  of 
the  amount  of  field  information  (RMS)  beyond  the  harmonic  degree  N: 


<5Agf(N) 


m 

Y~  Cnr 

n=N 


(5.3) 


The  sum  (5.3)  is  important  because  given  the  resolution  N  of  a  particular  field 
approximation,  the  sum  indicates  how  much  field  information  (RMS  anomaly)  is 
left  unresolved  by  the  approximation. 

A  third  method  of  global  characterization  is  obtained  by  giving  the  mean 
square  upward  continued  anomaly  on  a  spatial  sphere,  caused  by  boundary 
values  of  Ag  outside  a  cap  of  radius  i>ot  for  various  heights  of  the  spatial 
sphere.  This  mean  square  anomaly  is  the  same  as  the  mean  square  anomaly  of 
the  truncated  field  defined  by  (4.8)  and  can  be  obtained  from  (4.20)  by  using 
only  the  second  terra  of  that  formula  and  starting  the  summation  from  n=2: 


SAgrd,  fo)  -  q£(t,  fo)CnR 

n-2 

-IP) 


(5.4) 
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For  a  second  application  of  (4.20)  we  introduced  the  Rapp-180  field  to 
account  for  the  remote  zone  effects  (Nmax-180).  In  addition  we  used  a  more 
optimistic  set  of  errors  associated  with  the  Rapp-180  potential  coefficients, 
found  by  using  mis?=10  mgals  instead  of  20  mgals  in  equation  (30)  of  Rapp 
(1981).  Again  using  the  Tscherning/Rapp  Cn-model  and  upward  continuation 
distance  of  30  km,  the  total  error  (commission  plus  omission)  incurred  by  using 
the  Rapp-180  field  to  account  for  remote  zones  were  computed  by  (4.20)  for 
various  truncation  cap  radii  and  plotted  as  the  bottom  curve  of  Figure  2.  We 
see  that  now  the  total  error  is  drastically  reduced  compared  with  the  case  of 
the  top  curve  where  no  spherical  harmonic  model  is  used  to  account  for  remote 
zones.  For  a  truncation  radius  of  3*,  the  error  reduced  from  RMS  0.9  mgal  to 
RMS  0.1  mgal.  For  a  truncation  radius  of  1*,  the  bottom  curve  of  Figure  2 
shows  that  the  use  of  the  Rapp-180  field  to  account  for  the  remote  zone  field 
incurs  a  commission  plus  omission  error  of  0.45  mgal  global  RMS. 

Note  that  Figure  2  represents  a  global  error  analysis  that  may  not  be 
representative  of  a  local  area.  Figure  1  is  more  suited  for  analysis  of  specific 
areas.  For  example,  Figure  1  says  that  for  a  truncation  radius  of  3*  it  is 
immaterial  whether  we  use  the  Rapp-180  field  to  account  for  the  remote  zone  or 
not  because  the  Rapp-180  truncated  contributions  are  small  anyway.  On  the 
other  hand  Figure  2  says  that  not  using  the  Rapp-180  field  at  truncation  angle 
of  3*  causes  an  ommission  error  of  0.9  mgal  RMS,  and  the  use  of  Rapp-180  field 
decreases  this  error  to  0.1  mgal  (commission  plus  omission  error). 
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7.  Operational  Program  for  the  Generation  of  Anomaly  Field s 

In  this  section  we  summarize  the  programs  that  can  be  used  to 
operationally  generate  the  different  types  of  anomaly  fields  we  use  in  this 
report.  These  anomalies  are  the  Poisson  anomaly,  the  spherical  harmonic: 
anomaly,  the  topographic  anomaly,  and  the  Fourier  anomaly. 

What  we  call  the  Poisson  anomaly  field  is  generated  in  planar  approximation 
as  (2.3.1).  A  Poisson  anomaly  field  can  be  operationally  generated  using  a 
FORTRAN  program  fully  documented  in  Rapp  (1966).  The  program  accepts  5’x5* 
mean  anomalies  as  boundary  values  in  the  flat-earth  upward  continuation.  For 
a  more  detailed  field  generation  smaller  blocks  of  2.’5x2.'5  mean  anomalies  may  be 
input  near  the  computation  points,  and  the  program  then  automatically  rejects 
any  5’x5’  mean  anomalies  covered  by  the  input  2/5x2!5  values.  The  program 
also  has  the  useful  feature  of  rigorously  computing  the  normal  gravity 
corresponding  to  the  anomalies  being  computed,  the  sum  of  these  two  quantities 
being  a  model  for  observed  total  gravity.  The  program  exists  as  program  F499 
in  the  O^U  program  library. 

For  the  generation  of  spherical  harmonic  fields,  there  are  two  types  of 
existing  operational  programs  that  can  be  used.  If  the  intention  is  to  generate 
anomaly  values  at  individual  points  not  on  a  grid,  the  program  described  in 
Rapp  (1982)  can  be  used.  If,  however,  values  on  a  limited  grid  are  desired 
then  the  program  described  in  Rizos  (1979)  can  be  used.  The  program  by 
Rapp,  to  generate  values  at  individual  points,  exists  as  program  F477  in  the 
OSU  program  library,  and  the  program  by  Rizos,  to  generate  values  on  a  grid 
exists  as  program  F388.  A  comparison  of  these  and  other  spherical  harmonic 
programs  is  given  in  Tsc.herning  et  al.  (1983).  The  input  to  F477  and  F388  are 
the  set  of  potential  coefficients  and  the  geodetic  latitude,  longitude  and  height 
above  the  reference  ellipsoid  ($,  X,  h)  of  computation  points.  The  set-up  of  the 
programs  is  to  implement  equation  (2.2.2);  however,  for  reasons  stated  below 
equation  (2.2.2),  we  have  slightly  modified  F477  and  F388  for  our  applications 
to  compute  (2.2.3)  instead. 

Another  anomaly  field  f  interest  to  us  is  the  topographic  anomaly  field, 
generated  by  intergrating  the  gravitational  effects  of  topographic  masses  of 
assumed  density.  The  operational  generation  of  a  topographic  anomaly  field 
can  be  done  using  the  FORTRAN  program  described  in  Forsberg  (1984)  and 
existing  as  OSU  program  489.  There  are  various  modes  under  which  the 
program  runs,  as  detailed  in  Forsberg  (ibid.),  but  the  most  important  one  for 
the  purpose  of  our  studies  is  that  for  computing  the  external  gravity  anomaly 
field  generated  by  the  (positive  and  negative)  residual  masses  lying  between 
the  actual  topography  defined  by  a  digital  elevation  model,  and  a  reference 
topography  such  as  the  topography  to  180  spherical  harmonic  expansion. 
Another  mode  of  interest  to  the  procedures  recommended  in  this  report  is  the 
terrain  correction  (tc)  computations  which  will  be  needed  in  case  the  tc  are  not 
given  on  the  gravity  data  records. 

Finally,  by  Fourier  anomaly  we  mean  an  anomaly  which  is  upward  continued 
using  Fourier  transform  techniques.  For  the  generation  of  a  Fourier  anomaly 
field  we  used  a  simple  program  given  in  Appendix  B  and  existing  as  OSU 
program  F498.  The  theory  behind  this  program  is  detailed  in  Section  2.4  and 
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tests  are  conducted  in  Section  9.3.  A  more  extended  program  for  Fourier 
upward  continuation  exists  in  the  OSU  Department  of  Geology  (R.  von  Freest 
private  communication),  and  this  program  has  options  to  choose  from  among 
several  different  windowing  techniques.  However,  our  teeth  with  this  program 
have  not  shown  a  need  for  the  application  of  windowing  for  gravity  anomaly 
upward  continuation. 
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8.  Data  Preparation  for  Upward  Continuation  Teats  in  New  Mexico 

For  testing  alternative  models  for  the  operational  upward  continuation  of 
surface  free-air  anomalies,  we  prepared  data  in  a  7'x9*  area  in  New  Mexico. 
The  location  of  the  area  is  such  that  it  centers  the  site  of  the  balloon-borne 
gravity  project  coordinated  by  the  AFGL  (Lazarewicz,  et  al.,  1983),  with  surface 
data  coverage  extending  300  km  on  all  sides  of  the  balloon  flight. 


8.1  Available  Gravity  and  Elevation  Data 

Our  gravity  and  elevation  data  was  based  on  that  supplied  to  us  earlier 
(April,  1983)  by  the  National  Geodetic  Survey  (NGS).  The  gravity  data  were  in 
the  form  of  irregularly  distributed  point  values  of  surface  free-air  anomalies, 
as  shown  in  Figure  6.  The  gravity  data  record  included  the  following  items 
(Hittelman  et  al.,  1982): 

(1)  Geodetic  latitude  (♦),  geodetic  longitude  (X),  and  orthometric  height  (H)  of 
the  station. 

(2)  Measured  gravity  (g)  referenced  to  a  recoverable  base  station.  Base 
stations  had  been  adjusted  to  the  International  Gravity  Standardization 
Network  1971  (Morelli,  et  al.,  1972). 

(3)  Normal  gravity  (?)  at  the  Geodetic  Reference  System  1967  (GRS67)  reference 
ellipsoid,  computed  as: 

7  =  978031.85(1  +■  0.005278895  sin2  ♦  + 

+  0.000023462  sin4  ♦)  mgals  (8.1.1) 


(4)  Surface  free-air  anomaly  (Aga)  computed  as: 

Ags  =  g  +  0.3086  H  -  y  (8.1.2) 

where  0.3086  mgal/m  is  the  normal  gradient  of  gravity. 

(5)  Simple  Bouguer  anomaly  (Ag'n)  computed  as: 

Ag’e  =  Ags  -  0.1119  H  (8.1.3) 

where  the  term  0.1119  H  is  the  attraction  of  a  Bouguer  plate  of  standard 
continental  density  of  2.67  g/cm3. 

(6)  Terrain  correction  (tc)  formally  given  by  (Moritz,  1966,  p.  88): 
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Figure  6.  Point  Location  of  Original  NGS  Gravity  Data  Before  Thinout 
Procedure,  New  Mexico  Balloon  Gravity  Test  Site 


(8.1.4) 


where  k  =  Newtonian  gravitational  constant,  p  -  density;  R  =  mean  earth 
radius;  H  =  elevation;  P  =  point  to  which  tc  refers;  a  -  unit  sphere; 
da  =  element  of  solid  angle;  l0  :  2  R  sin  i>/2\  i*  -  angular  distance  between 
P  and  da. 

(7)  Standard  error  of  Ags  and  Ag’3, 

(8)  Various  codes:  Agency  code,  quality  code,  elevation  code,  and  source 

code. 


There  were  a  total  of  18,386  original  gravity  points  in  our  area.  Out  of  these, 
13,455  were  NGS  coded  as  ’ACCEPTED’  while  3,931  were  coded  as  'NOT  EDITED’. 
We  decided  to  consider  all  the  18,386  points  regardless  of  code  as  input  to  the 
data  thinning  step  (see  Section  8.2  below). 

The  elevation  data  were  in  the  form  of  30x30  arcsec  grid  point  values.  The 
data  covered  our  7*x9*  New  Mexico  area  except  at  three  l*xl*  blocks  in  the 
southwest  corner  from  latitude  29*  to  30*  and  longitude  251*  to  254*.  We 
decided  to  fill  these  missing  blocks  by  5’x5’  mean  elevations  from  data  supplied 
by  the  Defense  Mapping  Agency  (DMA).  The  NGS  data  were  used  in  our 
procedures  both  as  the  original  30"x30"  point  values  and  to  obtain  5’x5’ 
averaged  values.  The  5’x5’  mean  elevations  were  formed  by  straight  averaging 
all  30"x30"  values  that  fell  inside  the  5’x5’  block  (disjoint  averages).  A 
contour  map  of  the  topography  in  our  area  based  on  5'x5*  mean  elevations  is 
shown  in  Figure  7  for  a  contour  interval  of  50  meters. 


8.2  Data  Thinning 

The  original  set  of  18,386  NGS  points  shown  in  Figure  6  were  input  into  a 
data  thinning  procedure.  This  step  was  done  to  make  the  data  distribution 
more  uniform  and  to  later  avoid  collocation  inversion  problems  associated  with 
data  points  that  are  very  close  together.  To  thin  out  the  data  a  single  pass 
was  made  to  select  only  the  first  point  that  fell  inside  each  element  of  a 
3.5  km  3.5  km  (A*  =  2.’0,  AX  =  215)  grid  mesh.  After  the  thin  out  procedure  a 
total  of  10,208  data  points  were  left.  Of  these,  2  points  were  later  discovered 
as  blunders  and  removed  leaving  a  final  selection  of  10,206  points  which  are 
shown  in  Figure  12  of  Section  8.7. 


8,3 _ 5’x5’  Mean  Anomalies 

A  mean  anomaly  can  be  derived  by  first  subdividing  the  mean  block  into 
pxq  sub-blocks,  predicting  point  anomalies  at  the  centers  of  the  sub-blocks, 
then  averaging  the  predicted  point  anomalies.  A  predicted  point  anomaly  is 
given  by  least  squares  collocation  as  follows:  (Rapp,  1978,  p.  134): 


-  • 
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4Si  -  cij  (cjj  f  Djj)  1  4gj 


(8.3.1) 


where 


C^j  anomaly  covariance  vector  between  the  point  i  being  predicted 

and  the  data  points  j. 

Cjj  anomaly  covariance  matrix  of  the  data  points 

Dj j  error  covariance  matrix  of  the  data  points,  taken  as  diagonal 

with  elements  the  variances  of  the  data  points. 


Aspects 

known 

(8.6). 


on  the  covariance  function  to  use  in  (8.3.1)  and  on  the 
trends  in  Agj  are  discussed  separately  in  Sections  (8.4), 


removal  of 
(8.5),  and 


The  mean  anomaly  is  related  to  the  pxq  center  point  values  inside  the 
block  by 


l-l 


Substituting  (8.3.1)  into  (8.3.2)  we  get 


A  g  =  C 


Agj 


(cjj 


Agj 


where 


i  =  l 


(8.3.2) 


(8.3.3) 


(8.3.4) 


is  approximately  the  covariance  vector  between  the  mean  value  and  the  data 
points  j.  Equation  (8.3.4)  expresses  a  numerical  integration  procedure 
(Heiskanen  and  Moritz,  1967,  p.  277)  for  the  determination  of  mean  value  to 
point  value  covariance  using  a  i^cq  subdivision  of  the  mean  block.  Equation 
(8.3.3)  expresses  the  direct  prediction  of  the  mean  anomaly  from  given  point 
anomaly  data.  The  standard  error  of  the  predicted  is  the  square  root  of: 


m3~  r  C 
Ag 


AgAg 


-  C—  ;T  (Cjj  +  Djj) 


Agj 


Agj 


(8.3.5) 


where  C—  —  is  the  variance  of  the  mean  value  being  predicted,  given 
approximately  by 
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(8.3.6) 


r -  - - = — 

4«4g  (Pxq)2 

i=l  k=  1 


that  is,  is  the  average  of  all  covariances  between  subdivision  block 

center  points  inside  the  mean  block. 

For  our  applications  to  be  detailed  later  (Section  8.7)  we  used  equations 
(8.3.3)  to  (8.3.6)  to  predict  5’x5’  mean  anomalies  from  the  the  thinned  out  data 
of  Section  8.2.  We  used  a  2x2  block  subdivision  (p-2,  q=2)  and  only  the  ten 
closest  data  points  to  the  center  of  the  5’x5’  block  being  predicted.  The 
limitation  to  ten  data  points  was  motivated  by  the  computer  expense  required 
to  invert  the  matrix  in  (8.3.3)  with  dimension  equal  to  the  number  of  data 
points.  This  (approximate)  collocation  from  the  ten  closest  data  points  was 
compared  with  a  more  rigorous  but  much  slower  collocation,  giving  a  mean 
difference  of  0.1  mgal  and  an  RMS  difference  of  2.7  mgalB  for  a  12x12  array  of 
prediction  points  in  the  data  sparse  l*xl*  area  from  latitude  30*  to  31*  and 
longitude  254*  to  255*  (See  Figure  12  for  data  distribution).  The  predicted 
quantities  were  refined  Bouguer  anomalies  with  roughness  (standard  deviation 
from  the  mean  =  15  mgals  for  the  l'xl*  test  area)  shown  in  Figure  13.  The 
differences  between  the  rigorous  and  approximate  collocation  are  expected  to 
get  smaller  in  the  immediate  area  of  the  balloon  flight  because  of  the  increased 
density  of  data  there.  The  rigorous  collocation  used  a  one-time  inversion  of  a 
matrix  of  size  524x524,  giving  data  coverage  out  to  0:5  away  from  the 
prediction  points  while  the  approximate  collocation  used  144  inversions,  with 
each  inversion  involving  a  10x10  matrix.  The  rigorous  collocation  was  about 
100  times  slower  than  the  approximate  collocation  in  this  test  case. 

Theoretical  and  practical  aspects  on  the  choice  of  covariance  functions  to 
use  in  (8.3.3)  and  ways  to  de-trend  the  data  are  discussed  in  the  next  three 
sections. 


8.4  Covariance  Function 

The  covariance  function  to  use  in  (8.3.1)  or  (8.3.3)  is  well-defined  only  in 
the  global  case.  However,  for  local  applications  the  choice  of  covariance 
function  to  use  is  rather  arbitrary.  This  arbitrariness  comes  from  the  fact 
that  the  covariance  function  to  use  should  be  tailored  to  approximate  the  local 
empirical  covariance  function,  and  this  local  function  has  no  clear-cut 
theoretical  definition.  The  most  important  questions  are:  what  size  of  local 

area  should  be  used  to  derive  the  function?,  and  what  trends  should  be 

removed  from  the  data?  Once  a  practical  decision  has  been  made  on  these  two 

items,  however,  the  estimation  procedure  becomes  clear:  The  empirical 

covariance  function  is  first  derived  by  averaging  products  of  de-trended  data 
sumples  in  the  specified  area,  according  to  the  definition  of  covariance  function 
(Heiskanen  and  Moritz,  1967,  p.  253).  The  derived  empirical  covariance  function 
is  then  usually  approximated  by  a  simple  analytical  function  that  lends  itself  to 
closed  form  covariance  propagation  if  needed,  resulting  in  a  self-consistent 
system  for  the  estimation  of  various  linear  functionals  of  the  earth’s  anomalous 
potential.  The  way  to  approximate  an  empirical  covariance  function  by  an 
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analytical  function  of  specified  form,  is  by  satisfying  the  three  essential 
parameters  of  the  empirical  function,  namely,  the  variance,  correlation  length, 
and  curvature  parameter  (Moritz,  1980,  p.  174). 

For  our  applications  we  used  a  simple  tailoring  procedure  for  the  covar¬ 
iance  function.  We  took  the  global  anomaly  covariance  function  of  Tscherning 
and  Rapp  (1974)  and  first  subtracted  the  first  36  harmonics.  This  resulted  in 
a  new  covariance  function  with  a  correlation  length  of  about  20  km,  approx 
imating  the  correlation  length  of  the  de-trended  data  (see  below)  used  in  out- 
estimation  procedures.  The  covariance  function  was  then  scaled  to  satisfy  the 
variance  of  the  de-trended  data.  This  scaling  would  not  affect  the  previously 
tailored  correlation  length.  For  the  curvature  parameter,  there  was  no  specific 
treatment  given  to  satisfy  the  data;  the  practical  problems  associated  with  the 
empirical  computation  of  the  curvature  parameter,  as  well  as  its  approximate 
computation  by  finite  differences,  can  be  found  in  Schwarz  and  Lachapelle 
(1980). 


8 . 5_  Data  De-Trending  at  High  Frequencies 

It  is  well-known  that  short  wavelength  free-air  anomalies  are  strongly 
correlated  with  short-wavelength  topography  -  for  an  instructive  physical 
interpretation  of  this  fact  using  a  simple  crustal  density  and  isostatic  com 
pensation  model,  see  Moritz  1968,  p.  28.  Therefore,  the  computational  removal 
of  the  attraction  caused  by  topographic  masses  is  certain  to  remove  most  of 
the  roughness  that  may  be  present  in  an  anomaly  field. 

This  removal  of  roughness  of  the  field  is  very  important  in  interpolation 
problems.  With  the  removal  of  as  much  roughness  as  possible,  the  correlation 
length  of  the  residual  field  will  be  enlarged  as  much  as  possible.  This  implies 
that  the  ratio 

correlation  length 

P  mean  data  spacing 


will  also  be  enlarged  as  much  as  possible,  and  this  is  a  key  to  strengthening 
the  interpolation  of  the  residual  field  (see  Sunke),  1981,  pp.  88-93,  where 
values  of  p  of  at  least  p=3  are  indicated  to  be  desirable).  Of  course,  once  the 
interpolated  value  from  the  residual  field  has  been  obtained, the  total  field 
value  can  be  obtained  by  adding  back  the  influence  of  topographic  masses. 
This  influence  is  computable  from  detailed  topographic  height  data  which  are 
assumed  to  be  readily  available  in  the  prediction  area. 

Anomaly  data  de-trended  at  high  frequencies  can  be  the  simple  or  the 
refined  Bouguer  anomalies.  The  simple  Bouguer  anomaly  is  given  by  (8.1.3): 


i g ’ -  ag.  -  0.1119  H  (8.5.  I  ) 


i 


i 


4 
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where  -0.1119  II  means  the  removal  from  surface  data  Ags  of  the  gravitational 
attraction  of  a  moving  Bouguer  plate  of  standard  density  2.67  g/cm3.  The 
refined  Bouguer  anomaly  is  given  by: 

AgB  =  Ags  “  0.1119  H  +■  tc  (8.5.2) 

where  tc  is  the  gravimetric  terrain  correction  formally  given  by  (8.1.4). 
Whereas  (8.5.1)  represents  the  removal  of  a  moving  Bouguer  plate,  the 
application  of  tc  to  (8.5.1)  to  arrive  at  (8.5.2)  means  that  now  in  (8.5.2)  the 
gravitational  attraction  of  the  actual  (non-moving,  fixed)  topography  is  removed 
from  Ags.  The  tc  were  given  for  our  NGS  gravity  data;  if  they  had  not  been 
given,  we  would  have  had  to  compute  them  using  the  operational  program  by 
Forsberg  (see  Section  7).  Such  computations  of  tc  tend  to  be  expensive  (about 
0.2  CPU  sec  per  point  on  the  Amdahl  470  V/6  system)  if  they  have  to  be  done 
for  all  observation  points.  In  this  respect  the  studies  by  Sideris  (1984)  on  the 
computation  of  tc  by  Fast  Fourier  Transforms  (FFT)  should  prove  to  be  very 
important. 


8.6  Data  De-Trending  at  Low  Frequencies 

The  Bouguer  anomalies  produced  by  (8.5.1)  or  (8.5.2)  are  much  smoother 
than  the  original  Agg  field  but  are  biased,  having  large  and  systematically 
negative  values  in  mountainous  areas  -  again,  it  is  instructive  to  see  Moritz 
(1968,  p.  28)  for  a  physical  explanation  of  this  fact  using  isostatic  compensation 
theory.  In  accordance  with  the  statistical  aspect  of  the  least  squares 
collocation  interpolation  procedure,  gross  trends  should  first  be  removed  from 
the  data  before  interpolation  (see  Moritz,  1980,  Sec.  38:  "The  Meaning  of 

Statistics  in  Collocation").  To  de-trend  the  Bouguer  anomaly  data,  one  way 
would  be  to  postulate  a  low  order  trend  surface  and  fit  this  surface  to  the 
data,  possibly  in  the  context  of  least  squares  collocation  with  systematic 
parameters  (Sunkel,  1983),  or  possibly  in  the  context  of  a  simple  least  squares 
adjustment.  For  our  purposes  we  de-trended  the  Bouguer  anomaly  data  using 
available  spherical  harmonic  expansions  to  180  of  free-air  anomaly  and  topo¬ 
graphy.  The  effect  that  the  de-trending  at  low  frequencies  has  on  inter¬ 
polation  results  is  described  later  under  Section  8.7,  step  (7). 

The  free-air  anomaly  on  an  equatorial  sphere  can  be  generated  from 
potential  coefficients  to  degree  Nmax,  using  (2.2.3)  with  Ho=0: 

,  \ax  n  _  _  _ 

Ags  =  ^  (n-1)  ^  (C^n,  cos  mX  +  Snn)  sin  raX)  Pnm  (sin  <t>)  (8.6.1) 

n  -2  ra-0 


(the  superscript  s  denotes  spherical  harmonics).  The  topography  can  also  be 
expanded  in  terms  of  spherical  harmonics  to  degree 


Hs 


^max 

E 

nr0 


E 


m=0 


( Afj,,,  cos  mX  +  Bnra  sin  niX)  Pnm 


(sin  ♦) 


(8.6.2) 


We  then  have  a  Bouguer  anomaly  to  resolution  corresponding  to  harmonic 
degree  Nmax: 


Ag§  r  AgS  _  0.1119  Hs 
B 


(8.6.3) 


Subtraction  of  the  reference  Bouguer  anomaly  Agsg  from  the  Bouguer  anomaly 
in  (8.5.1)  and  (8.5.2)  gives,  respectively,  the  residual  Bouguer  anomaly  (8.6.4) 
and  terrain  corrected  residual  refined  Bouguer  anomaly  (8.6.5)  (see  equation 


(3.3.13)  for  a  fuller  understanding  of  equation  (8.6.5)): 

Agf’  =  (Ags  "  0.1119  H)  -  ( Ags  -  0.1119  Hs )  (8.6.4) 

(agr  r  tcs)  --  (Ags  -  0.1119  H  +  tc)  -  (4gs  -  0.1119  Hs)  (8.6.5) 

The  last  two  equations  can  also  be  written  as 

4gr’  ’  Ags  -  4gs  -  0.1119  (H  -  Hs)  (8.6.6) 

(4gr  *-  tcs)  -  Ags  -  Ags  -  0.1119  (H  -  HS)  +  tc  (8.6.7) 


The  last  two  equations  state  that  the  original  anomalies  Ags  are  de-trended  (1) 
in  the  long  wavelength,  by  subtracting  free-air  anomalies  Ags  generated  from 
spherical  harmonic  expansion;  and  (2)  in  the  short  wavelength,  by  doing 
"Bouguer  reduction"  not  with  respect  to  the  geoid  but  with  respect  to  the 
higher  order  but  still  smooth  surface  Hs  from  spherical  harmonics. 


8.7  Actual  Predictions  of  5’x5’  Mean  Anomalies 


In  accordance  with  the  previous  discussions  we  took  the  following  steps  to 
predict  5’x5'  mean  anomalies,  for  use  in  our  anomaly  upward  continuation 
procedures.  The  starting  anomaly  data  were  the  point  surface  free-air 
anomalies  (look  ahead  to  Figure  12  for  point  location),  resulting  from  the  thin 
out  procedure  of  Section  8.2.  As  stated  in  Section  8.1  the  gravity  record  also 
contained  the  elevation  H  of  the  station  and  the  terrain  correction  tc.  Hero  are 
the  various  steps: 

(1)  Referejice  free-air  anomalies  Ags  (8.6.1)  were  generated  on  a  01 2 5x0". 2 5  grid 
using  the  Rapp- 180  (Rapp,  1981)  potential  coefficients  with  Nlmnx=]80.  A  con 
tour  map  of  this  data  is  shown  in  Figure  8  with  a  contour  intervul  of  5  mgals. 
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Figure  8 


Rapp- 180  Gravity  Anomalies  at  Zero  Altitude,  New  Mexico 
Gravity  Test  Site.  C.I.  =  5  mgals. 


Reference  elevations  Hs  (8.6.2)  were  generated  on  a  0".  2  5x0".  2  5  grid  using 
iographic  coefficients  available  at  OSU  (tape  GS140,  file  15)  with  Nmnx=lHO.  A 
itour  map  of  this  data  is  shown  in  Figure  9  with  a  contour  interval  of  50 
ters. 

Reference  Bouguer  anomaly  values  Agsn  (8.6.3)  were  generated  on  a 
5x0*.25  grid  from  the  Ags  and  Hs  of  steps  (1)  and  (2).  This  data  is  contoured 
Figure  10  with  a  contour  interval  of  5  mgals.  In  our  procedures  we  could 
o  generate  AgSg  directly  in  one  step  because  we  had  combined  the  two 
■ies  in  (8.6.1)  and  (8.6.2)  to  produce  a  Bouguer  anomaly  series  to  degree  180, 
,h  its  own  Bouguer  anomaly  spherical  harmonic  coefficients. 

Refined  Bouguer  anomalies  Agn  (8.5.2)  were  computed  at  the  irregularly 
itributed  data  points  using  A gs,  H  and  tc  given  on  the  gravity  records.  The 
S  value  of  the  original  irregularly  distributed  Ags  was  26  mgals.  The  RMS 
lue  of  the  refined  Bouguer  anomalies  Agp  with  the  mean  removed  was  higher, 
mgals,  because  although  Agg  was  smooth  it  had  significant  long  wave  trend. 

Terrain  corrected  residual  refined  Bouguer  anomalies  Agr  +  tcs  (8.6.5)  were 
mputed  at  the  irregular  data  points  by  first  interpolating  the  0t25xCr“.25  grid 
Agsg  from  step  (3)  to  obtain  the  reference  Bouguer  anomaly  at  the  data 
int,  then  subtracting  this  reference  value  from  the  refined  Bouguer  anomaly 
step  (4). The  RMS  value  of  the  irregularly  distributed  anomalies  (Agr  +  tcs) 
is  a  smooth  and  centered  15  mgal. 

)  5’x5’  mean  values  of  terrain  corrected  residual  refined  Bouguer  anomalies 
.6.5)  were  predicted  from  the  data  of  step  (5)  using  the  "collocation  from  the 
jsest  10  points"  procedure  described  in  Section  8.3.  To  repeat,  a  2x2 
bdivision  of  the  5’x5’  block  was  used.  The  covariance  tailoring  procedure 
;ed  is  described  in  Section  8.4.  A  contour  map  of  the  predicted  5’x5’  mean 
sidual  refined  Bouguer  anomalies  is  shown  in  Figure  11,  with  a  contour 
terval  of  5  mgals.  Note  that  this  de-trended  anomaly  surface  is  much 
mother  and,  therefore,  much  more  reasonable  to  interpolate  than  the  trended 
iginal  Ags  surface  (look  ahead  to  Figure  14).  The  point  location  of  the 
regularly  distributed  data  from  which  5’x5’  predictions  were  made,  as  well  as 
e  formal  standard  errors  of  predicion  coming  out  of  the  collocation 
■ocedure,  are  shown  in  Figure  12.  The  collocation  error  estimates  shown  are 
least  correct  on  a  relative  basis,  but  also  their  absolute  values  are  expected 
be  meaningful  because  of  the  use  of  an  empirically  tailored  covariance 
notion  (Schwarz  and  Lachapelle,  ’980,  p.  33). 

)  A  "back  solution"  could  now  be  made  starting  from  the  predicted  5’x5’ 
?an  values  of  step  (6).  5’x5’  grid  point  values  of  AgsB  (8.6.3)  were  first 

terpolated  at  the  centers  of  the  5’x5‘  mean  blocks  of  step  (6)  using  the 
25xCT.25  grid  of  Agsg  values  from  step  (3).  These  interpolated  values  were; 
en  added  to  the  values  in  (8.6.5)  from  step  (6)  to  produce  fFx5_’  meanvalues 
refined  Bouguer  anomalies  Agg  (8.5.2),  contoured  in  Figure  13  with  a  contour 
terval  of  5  mgals.  The  difference  between  the  Bouguer  anomalies  in  Figures 
I  (from  this  step)  and  11  (from  step  (6))  is  the  presence  of  long  wave 
luguer  anomaly  trend  in  Figure  13. 

To  test  the  significance  of  de-trending  at  low  frequencies  we  repeated  the 
■edirtions  of  the  quantity  shown  in  Figure  13,  this  time  using  the  simple  data 
■  erage  as  reference  value,  instead  of  the  "Bouguer  anomaly  to  180".  The 
(variance  function  of  step  (6)  was  accordingly  scaled  to  reflect  the  variance 
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le  6.  Statistics  of  Upward  Continued  Terrain  Correction,  (tc,)[{0 


Altitude 
of  Upward 
Continuation  (H0) 

Mean 

Std.  Dev. 

RMS 

28.5  km 

0.38  mgal 

*0.61  mgal 

0.71  mgal 

8.5 

0.52 

*0.65 

0.82 

3.5 

0.70 

±0.66 

0.94 

0.0 

0.72 

±1.16 

1.32 

.■  relatively  fast  decrease  of  standard  deviation  from  *1.16  mgals  (Hoz0  km) 
*0.66  mgals  (Ho^.S  km),  illustrates  that  tc,  has  energy  in  the  very  high 
quency  range,  and  this  energy  gets  lost  by  attennuation  at  a  very  short 
ward  continuation  distance  H0. 

As  stated  in  Section  8.7,  the  terrain  correction  has  a  very  Bhort  correlation 
gth,  but  it  also  contains  a  weak  long  wavelength  signal.  This  long 
velength  signal  attenuates  rather  slowly  with  altitude,  and  shows  up  in  Table 
is  some  significant  effect  (RMS  0.71  mgal)  even  at  28.5  km  altitude. 

.2  Direct  vs.  Indirect  Upward  Continuation  Terms 

To  analyze  further  the  numerical  differences  between  the  direct  and 
iirect  upward  continuation  methods,  let  us  first  review  what  we  did  in  the 
lirect  method.  First  we  performed  the  following  split  of  the  surface  free-air 
:>malies: 

4gs  --  4gr  ,  4gs  f  4gt  (9.2.4) 


ere  Ag8  is  the  spherical  harmonic  contribution,  ig*-  is  the  influence 
ntributed  by  topographic  masses,  and  Agr  is  the  residual  anomaly.  Then  to 
th  sides  of  (9.2.4)  we  added  the  terrain  correction  tcs  of  the  reference 
xigraphy  HH  (for  Hs  we  used  a  topography  to  spherical  harmonic  expansion 
1): 

figs  *•  tcs)  Ugr  t-  tc3)  4gs  +  Agt  .  (9.2.5) 
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ure  15a.  Terrestrial  Anomaly  Profile  and  its  Uplifted  Version  at  28.5  km 
Altitude,  New  Mexico  Balloon  Gravity  Test  Site. 

Profile  Latitude  =  32".5417. 


ure  15b.  Rapp-180  Anomaly  Profile  and  its  Uplifted  Version  at  28.5  km 
Altitude,  New  Mexico  Balloon  Gravity  Test  Sito. 

Profile  Latitude  =  32*.5417. 
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into  the  causes  of  the  numerical  differences  between  the  direct  and  indirect 
methods. 

For  visualization  Figure  15a  shows  the  60-point  terrestrial  anomaly  profile 
and  its  upward  continued  version  at  altitude  H0=28.5  km.  Correspondingly  the 
Rapp- 180  anomaly  profile  at  the  ground  level  and  at  altitude  28.5  km  are  shown 
in  Figure  15b. 


9.2  Other  Studies 


This  section  presents  other  studies  that  we  conducted  using  the  data  in 
New  Mexico.  The  objective  is  to  give  more  information  on  the  numerical  aspects 
of  various  procedures  related  to  upward  continuation  of  surface  free-air 
anomalies. 


9.2.1  Formally  Upward  Continued  Terrain  Correction 

For  the  direct  method  we  had  the  ter  rain -uncorrected  version: 

(ags)Ho  ^  Up{4gs}  =  ^  dxdy  (9.2.1) 

A  0 

and  the  terrain-corrected  version: 

(4gs  +  tc)j}0  =  Upfags  +  tc}  =  II  p3tC^  tody  (9.2.2) 

A  D° 

The  numerical  difference  between  these  two  versions  of  the  direct  method  is: 
(4f?s  f  tcjHo  “  (4«s)Ho  -  ( tci  )8o  >  (9.2.3) 


that  iB,  the  upward  continued  terrain  correction.  The  subscript  "1”  is  used  in 
(9.2.3)  to  indicate  that  the  terrain  correction  tcj  is  not  equal  to  the  true 
terrain  correction  but  is  rather  a  terrain  correction  that  contains  errors 
caused  by  the  errors  in  the  individual  predicted  quantities  4gs  and  (4gs  +  tc). 
In  other  words  tc,,  formed  as  the  differences  between  predicted  4gK  and 
(4g8  t  tc),  can  be  called  a  "predicted"  terrain  correction.  A  comparison  of  tc, , 
with  the  rigorously  computed  tc  has  already  been  given  in  Section  8.7.  For 
the  upward  continued  tc,,  i.e.  (tc,)DHo  we  found  the  statistics  shown  in  Table 
6.  The  statistics  were  found  for  various  upward  continuation  distances  H0, 
using  the  60-point  upward  continuation  profile  described  in  the  previous 
section  and  5’x5'  tc,  data  in  the  7*x9*  New  Mexico  test  area. 
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Table  3.  Comparison  of  Upward  Continued  Anomalies.  Upward  Continuation  Distance  H„  -  28.5  kin. 
(Superscript  1:  Indirect  Method;  D:  Direct  Method). 


9.  Numerical  Investigations 


9.1  Comparison  of  Direct  and  Indirect  Upward  Continuation  Results 

We  applied  the  indirect  and  direct  upward  continuation  methods  described 
in  Section  3.3  to  our  7*x9*  study  area  in  New  Mexico.  Figure  7  on  page  52 
shows  the  topographic  features  in  the  area.  Upward  continuations  were  done 
for  a  60-point  profile  running  east-west  from  longitude  253*  to  258*.  Upward 
continuation  distances  used  were  H0=28.5,  8.5,  and  3.5  km.  These  rather  odd 
values  of  H0  resulted  from  considering  the  data  to  be  at  a  mean  elevation  of 
1.5  km  and  the  uplifted  profiles  to  be  at  30,  10  and  5  km  elevation.  The  60 
upward  continuatuation  points  were  located  at  the  centers  of  the  5’x5’  mean 
data  blocks  directly  beneath  the  profile.  Although  60  points  were  used,  it  was 
sufficient  to  present  results  only  for  every  6^  point  giving  11  presentation 
points. 

The  entire  84x108  grid  of  5’x5’  mean  anomaly  data  (Section  8.7)  and  5’x5’ 
mean  elevation  data  (Section  8.1)  covering  the  7*x9*  area  Were  used  in  upward 
continuations.  The  types  of  anomalies  used  were  the  terrain  corrected  residual 
refined  Bouguer  anomalies  (Agr+tc8)  for  the  indirect  method,  and  the  surface 
anomalies  Ags  and  terrain-corrected  surface  anomalies  (Ag8  +•  tc)  for  the  direct 
method.  The  detailed  1km  x  1km  point  elevation  data  (Section  8.1)  were  used 
near  the  computation  points  for  the  prism  integration  of  the  attraction  of 
topographic  masses  needed  in  the  indirect  method.  Also  for  the  purpose  of  the 
indirect  method,  the  Rapp-180  potential  coefficients  and  a  set  of  degree  180 
spherical  harmonic  coefficients  for  the  topography  (Section  8,6)  were  used  to 
generate  reference  values  of  gravity  anomalies  and  topography.  The  opera¬ 
tional  programs  used  in  the  numerical  investigations  are  described  in 
Section  7. 

Tables  3  to  5  give  the  results  of  the  indirect  and  direct  upward  continu¬ 
ations  for  the  11  presentation  points,  for  H0=28.5,  8.5,  and  3.5  km.  Columns  1 
to  3  of  the  tables  give  the  three  components  of  the  indirect  method,  namely  (1) 
anomaly  contribution  from  the  medium  wavlength  part  of  terrestrial  data,  (2) 
anomaly  contribution  from  the  long  wavelength  spherical  harmonic  field,  and  (3) 
anomaly  contribution  from  shallow  topographic  masses.  Column  4  gives  the 
total  anomaly  (Ags  +  tc8)  from  the  indirect  method  (see  page  20  for  the 
rationale  behind  the  use  of  tc8). 

The  results  from  the  indirect  method  were  compared  with  those  from  the 
direct  method,  with  differences  shown  in  columns  5  and  6  of  Tables  3  to  5. 
Column  5  shows  thut  the  direct  method  using  the  surface  anomalies  alone  (Ags) 
produces  a  profile  that  is  systematically  too  low  compared  with  the  expectedly 
more  rigorous  profile  of  the  indirect  method.  The  biases  have  mean  values  of 
0.64  mgal  (H0=28.5  km),  0.54  mgal  (H0=8.5  km),  and  0.71  mgal  (H0=3.5  km). 
Column  6  shows  that  the  use  of  terrain-corrected  surface  anomalies  (AgH4-tc)  in 
the  direct  method  improves  the  bias  of  column  5  down  to  0.03  mgal  (Ho -28.5 
km),  0.11  mgal  (H0=8.5  km),  and  0.06  mgal  (H0=3.5  km).  The  standard  deviation 
of  the  differences  of  columns  5  and  6  (direct  method)  with  column  4  (indirect 
method)  also  improves  slightly  with  the  use  of  (Ags+tc)  instead  of  Agg  in  the 
direct  upward  continuations.  In  the  next  section  we  will  make  further  studies 
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found  by  Forsberg,  1984,  p.  82).  A  theoretically  better  alternative  would  have 
been  to  take  the  Faye  anomalies  of  step  (8)  and  subtract  rigorously  computed 
tc  values  using  Forsberg’s  (Section  7)  program,  but  this  alternative  would  be 
prohibitively  expensive  at  the  present  time.  A  possible  way  out  of  this 
expense  would  be  the  development  of  FFT  techniques  to  compute  tc  (Sideris, 
1984). 

The  "predicted"  tc  of  step  (9)  was  brought  to  light  as  follows: 


tc( predicted)  (4gs  +  tc)STEP(8)  AgsSTEP(9) 


Values  of  tc(predicted)  were  taken  in  the  l*xl*  area  from  latitude  32*  and  33* 
and  longitude  254"  to  255*,  and  these  values  compared  with  rigorously 
computed  tc  from  Forsberg’s  program.  We  found  the  following  statistics  using 
a  6x6  grid  of  comparison  points  (units:  mgals): 


Statistics 


Difference 

tc(predicted)  -  tc( rigorous) 


Actual  Value 
tc( rigorous) 


Mean 

-1.18 

Std.  Dev. 

*2.69 

RMS 

2.06 

Maximum  Absolute  Value 

6.06 

2.68 

*2.95 

3.99 

14.38 


On  average,  one  can  say  from  the  above  table  that  50%  of  the  true 
tc-inforraation  has  been  recovered  in  the  prediction.  Inspection  of  the  actual 
differences  which  are  not  given  here  reveals  that  the  general  shape  of  the 
tc-function  can  be  reasonably  predicted,  but  that  a  reasonable  prediction  of 
detailed  features  basically  relies  on  the  chance  that  there  is  a  data  point  close 
to  the  prediction  point.  One  could  expect  that  for  a  reasonable  prediction  of 
detailed  features  of  tc  the  data  spacing  would  have  to  be  much  less  than  the 
correlation  length  of  tc  in  the  area,  say,  one-third  the  correlation  length;  this 
is  indicated  from  the  discussions  of  Sunkel  (1981,  pp.  88-93)  who  gives  general 
data  density  requirements  as  a  function  of  the  correlation  length  of  the 
function  being  interpolated. 
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of  the  "mean-value  centered"  data,  which  was  2400  mgal2.  For  comparison,  let 
us  distinguish  two  prediction  procedures  as  follows: 

Method  A:  Bouguer  anomaly  to  180  used  as  a  reference  surface.  Variance  of 
data  used  in  predictions:  225  mgal2. 

Method  B:  Simple  average  used  as  a  reference  value.  Variance  of  data  used  in 
predictions:  2400  mgal2. 

We  observed  the  following: 

(a)  In  areas  with  good  data  coverage  (i.e.,  areas  with  standard  error  of 
prediction  less  than  or  equal  to  5  mgals  as  shown  in  Figure  12),  the 
predicted  values  from  Methods  A  and  B  generally  agreed  to  better  than  0.2 
mgal. 

(b)  In  areas  with  poor  data  coverage  (standard  error  greater  than  or 
equal  to  10  in  Figure  12),  differences  of  7  mgals  were  observed  between 
Methods  A  and  B.  In  these  areas  the  used  reference  value  "anchors"  the 
predicted  value,  in  the  sense  that  the  predicted  value  tends  to  approach 
the  inference  value. 

(c)  The  most  significant  difference  between  Methods  A  and  B,  affecting 
both  areas  of  good  and  of  poor  data  coverage,  was  the  scaling  of  the 
formal  error  estimates.  Error  estimates  from  Method  B  were  more 
pessimistic  than  Method  A,  by  the  ratio  (2400/225)^. 

From  the  above  we  conclude  that  the  advantage  in  using  a  higher  order  trend 
surface  than  the  simple-average  surface  lies  in  a  better  scaling  of  the  error 
estimates;  the  predicted  values  themselves  are  not  critically  affected  by  the 
choice  of  trend  surface  for  areas  of  reasonable  data  coverage. 

(8)  To  the  5’x5’  refined  Bouguer  anomalies  (8.5.2)  from  step  (7)  we  added 
0.1119  H,  where  the  elevations  H  were  the  5’x5’  mean  elevations  mentioned  at 
the  end  of  Section  8.1.  In  accordance  with  (8.5.2),  the  results  would  be  5'x5’ 
mean  values  of  terrain  corrected  (i.e.  Faye)  surface  anomalies  (Ag8  +  tc), 
contoured  in  Figure  14  with  a  contour  interval  of  5  mgals.  Throughout  our 
procedures  we  of  course  assumed  that  5'x5’  mean  values  of  free-air  anomalies 
referred  to  5’x5’  mean  elevations,  and  this  assumption  is  justified  because  of 
the  strong  local  correlation  between  point  free-air  anomalies  and  elevations;  for 
a  further  discussion  of  correspondence  between  mean  free-air  anomalies  and 
mean  elevations,  see  Sunkel,  1981,  p.  5. 

(9)  It  was  also  of  interest  to  our  studies  to  predict  5’x5’  mean  anomalies 
without  huving  first  applied  the  terrain  corrections  to  the  irregular  data  Ags 
in  step  (4).  In  other  words,  we  essentially  repeated  all  the  de-trending, 
prediction,  and  back  solution  steps  working  with  simple  Bouguer  anomalies 
given  by  (8.5.2).  The  end  results  analogous  to  those  of  step  (8)  were  5'xS’ 
mean  yahies  of  _£ terrain,  uncorrected)  surface  free-air  anomalies  Agg. 

REMARK:  Note  that  implicitly  in  step  (9)  we  carried  the  quantity  "-tc"  (present 
in  the  original  Ags  data)  through  the  prediction  procedure.  Therefore, 
implicitly  we  attempted  to  predict  5’x5’  mean  values  of  "-tc”  which  would  have 
to  be  very  inaccurate  because  of  the  very  short  correlation  length  associated 
with  the  terrain  correction  (correlation  lengths  on  the  order  of  only  2.’5  were 
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Figure  13.  Refined  Bouguer  Anomalies,  New  Mexico  Balloon  Gravity  Test  Site 
C.I.  =  5  ragals. 
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Then  we  indirectly  upward  continued  (Ags  +  tcs) 


follows: 


a' 


(4gs  +-  tcs)Ho  =  (4£r  +  tcs)Ho  +  4SHo  +  4SHo 


(9.2.6) 


that  is,  the  indirectly  upward  continued  value  is  the  sum  of  three  separately 
upward  continued  terms:  the  first  term  is  a  direct  upward  continuation  by  the 
Poisson  integral;  the  second  term  is  an  upward  continuation  in  the  spherical 
frequency  domain;  and  the  third  term  is  an  upward  continuation  by  the  prism 
integration  of  the  gravitational  influence  of  topographic  masses. 

In  the  comparison  stage  (Tables  3,  4,  and  5)  we  get  a  good  mean  value 
agreement  between  the  set  of  results  from  the  indirect  method  and  those  from 
the  direct  method  that  used  the  terrain-corrected  surface  anomalies  (Ags  +  tc). 
The  latter  set  of  results  could  be  conceptually  obtained  by  adding  tc  to  both 
sides  of  (9.2.4)  then  applying  direct  upward  continuation: 


(4gs  +  tc)Ho  =  (4gr  +  tc)Ho  +  ( 4&s)Ho  +  (ASt)Ho 


(9.2.7) 


Therefore,  to  explain  the  numerical  differences  shown  in  Tables  3  to  5  between 
(Ags  +  tca)IHo  and  (iga  +  tc)DHo  we  need  to  examine  the  differences  between 
the  corresponding  terms  .  on  the  right  hand  side  of  equations  (9.2.6)  and 
(9.2.7). 

A.  Comparison  of  First  Terms 

The  difference  Aj  between  the  two  terms  is: 

A I  =  Ugr  +  tc)8o  -  (48r  +  tcS)g0  =  (tc  -  tcs)go  (9.2.8) 


As  stated  in  (9.2.3)  we  were  able  to  obtain  some  sort  of  predicted  tc  denoted 
by  tct.  On  the  other  hand,  the  determination  of  the  quantity  of  tc8  was 
problematic  in  terms  of  computer  time  requirements;  the  set-up  of  Forsberg’s 
program  that  we  were  using  (Section  7)  required  0.2  cpu  sec.  per  point  on  our 
AMDAHL  470  V/6  system,  and  we  needed  at  least  84x108-9072  values  of  tcs. 
Therefore,  we  did  not  specifically  perform  an  evaluation  of  Aj  =  (tc-tc8)^,,. 
However,  the  numerical  comparison  of  the  third  terms  below  show  that  on 
average  Aj  is  canceled  out  by  the  difference  (denoted  by  Ajjj)  between  the 
third  terms  (see  equation  (9.2.20)).  In  any  case,  we  would  expect  (tc-tc8)^H<> 
to  be  small,  because  the  long  wavelength  quantity  tc8  will  tend  to  cancel  the 
long  wavelengths  of  tc  (tc-tc8);  aF  we  have  seen  in  Table  6  only  these  long 
wavelengths  (and  not  the  very  short  ones)  would  have  had  the  chance  to  filter 
through  the  upward  continued  value  (tc-tce)^n0- 
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B.  Comparison  of  Second  Terms 


The  second  term  of  (9.2.7),  which  is  (4gs)^Ho>  uses  the  spherical  harmonic 
series  derived  Ags  and  upward  continues  this  using  the  Poisson  integral.  The 
second  term  in  (9.2.6)  conceptually  uses  the  same  quantity  Ags  but  upward 
continues  it  using  the  spherical  harmonic  series,  upward  continuation  being 
effected  by  multiplying  the  spectrum  of  Ags  by  the  factor  (R/r)n+2.  To 
compare  these  terms  we  look  the  Rapp-180  field  used  in  Section  9.1  and,  with 
program  F388  of  Section  7,  generated  5’x5’  center-point  values  of  Ags  on  the 
equatorial  sphere  level.  These  values,  taken  as  5’x5’  mean  values,  covered  the 
7*x9*  test  area  in  New  Mexico.  We  then  upward  continued  the  5’x5’  mean 
vnlues  by  Poisson  integration  (program  F499,  Section  7): 


“<•*.  -  t  i  I  f 

uo 

dxdy 

(9.2.9) 

and  compared  the  results 
(program  F388,  Section  7): 

with  the 

rigorously  upward 

continued  values 

o  D  kM  <r -  ,  .  . 

(*gs>Ho  =  2_  (n~l) 

f  a_ln+J 
la+H0J 

n 

^  (C^m  cos  mX  + 

n  =  2 

m=0 

Snm  sin  mX  )  Pnm  (sin  ♦)  (9.2. 10 ) 


The  implementation  of  equation  (9.2.9)  was  done  in  two  versions.  In  the 
first  version  (called  center-point  kernel  version)  the  contribution  of  a  5’x5’  Ags 
to  the  upward  continued  value  was  obtained  by  multiplying  Ags  by  the  step 
function  evaluation  of  the  integral  kernel  at  the  center  point  of  the  5’x5’  block. 
In  the  second  version  (called  integrated  kernel  version)  the  5’x5’  4  gs  was 
multiplied  by  the  rigorously  integrated  value  of  the  kernel  inside  the  area 

covered  by  the  5’x5’  block  (see  Rapp,  1966,  for  specific  equations).  The 

statistics  of  the  differences  a j I  =  (Ags)Dj.j0  “  iS'sHo  are  shown  in  Table  7.  The 

statistics  were  obtained  for  the  same  eleven  presentation  points  of  the  profile 

used  in  Section  9.1.  We  observe  that  the  integrated  kernel  implementation  of 
equation  (9.2.9)  can  keep  the  error  of  ( Ajgs ) DHo  under  1%  of  true  value 
This  is  true  down  to  the  lowest  upward  continuation  altitude  tested  which  was 
H"-l  km.  The  center-point  kernel  implementation,  on  the  other  hand,  can  keep 
the  relative  error  under  1%  down  to  10  km,  but  below  10  km  the  use  of  the 
integrated  form  is  necessary  since  the  tests  show  that  the  errors  of  the 
renter-point  form  reach  15%  at  5  km  and  blow  up  at  1  km. 

To  conclude,  the  second  term  AgsHo  of  equation  9.2.6  and  the  second  t<  rrn 
(4gH)U||,,  of  equation  9.2.7  are  in  close  agreement,  and  therefore  t'eir 

difference  4 [ [  is  not  a  major  contributing  factor  to  the  differences  found  in 
Section  9.1  (Tables  3  to  5)  between  the  direct  method  and  the  indirect  method 
of  upward  continuation.  However,  in  spite  of  the  small  differences  found 


between  Ag8Ho  and  (Ags)^Ho>  the  rigorous  quantity  igsHo  is  still  advisable  to 
use  (in  an  indirect  method  setting)  because  its  evaluation  does  not  increase  the 
computational  burden  very  much  and  it  has  the  advantage  that  in  principle  it 
uses  continuous  and  global  data  whereas  (Ag8)^Ho  (which  is  the  implicit 
evaluation  in  the  direct  method)  uses  a  step  function  approximation  of  the 
terrestrial  data  and  a  limitation  of  integration  to  within  a  finite  data  cap. 


Table  7.  Statistics  of  the  Absolute  Error  Incurred  in  Using  the  Poisson 
Integral  for  the  Upward  Continuation  of  the  Rapp-180  Anomaly  Field. 
5’x5’  Center  Point  Rapp-180  Anomalies  Used  as  Data,  Covering  the 
7*x9*  New  Mexico  Test  Area.  Units:  mgals. 


Altitude  of 

True  AgS 

Error 

Error 

Upward 

=Agjao 

(Center  Point 

( Integrated 

Continuation 

Kernel  Used) 

Kernel  Used) 

Ho=30  km 

Mean=  7.38 

0.05 

0.03 

S.D.=*6. 17 

*0.07 

RMS  =  9.43 

0.09 

0.08 

10  km 

9.02 

0.05 

*9.83 

*0.06 

*0.13 

13.01 

0.08 

0.15 

5  km 

9.52 

1.31 

0.10 

*11.05 

*1.63 

*0.13 

14.20 

2.03 

0.16 

1  km 

9.94 

107.46 

0.03 

*12. 15 

*132.98 

*0.03 

15.27 

166.20 

0.04 

C.  Comparison  of  Third  TermB 

Finally  let  us  now  turn  to  a  comparison  of  the  third  terms,  i.e.  ig^Ho  in 
(9.2.6)  and  (4gl)DHo  in  (9.2.7).  The  quantity  ag^Ho  is  the  indirect,  and 
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(AgM^Ho  the  direct,  upward  continued  value  of  Ag1.  jn  turn  4gt  j8  given  by 
equation  (8.7): 


Agt  =  0.1119  (H  -  Hs)  -  (tc  -  tc8) 


(9.2.11) 


where  in  our  case  the  reference  elevations  Hs  and  terrain  corrections  tcs  refer 
to  a  180-expansion  of  the  topography.  We  have  the  indirect  quantity: 


found  by  integration  of  the  attractions  caused  at  the 
computation  point  by  residual  topographic  masses  sub¬ 
divided  into  prism  integration  elements  (see 
Section  7  for  operational  program  used) ; 


and  the  direct  quantity: 


(a^Ho 


Ha  j  j  0.1119  (H  -  HS)  -  (tc  -  tc8)  ^ 
A  °o 


(9.2.13) 


As  we  indicated  under  equation  (9.2.8)  we  did  not  evaluate  tc8  because  of 
computer  time  limitations.  Therefore,  we  also  did  not  evaluate  (9.2.13)  where 
tc8  appears.  However,  tests  of  interest  could  still  be  conducted  by  dropping 
some  terms  in  (9.2.13);  we  used  two  versions: 

(Ag^Ho  -  2^  II  d*dy  (9.2.14a) 

A  D° 

with  Agt  =  Agt  -  tcs  =  0.1119  (H  -  Hs)  -  tc;  (9.2.14b) 


and 


(  Ag j  ^ {jo 


Ho 
2  n 


dxdy 


(9.2.15a) 


with  Agt 


t  - 


Ag1 


(tc  -  tcs)  =  0.1119  (H  -  HS) 


(9.2.15b) 


Using  our  5’x5’  data  in  New  Mexico,  various  upward  continuation  distances 
(H0=28.5,  8.5,  3.5  km),  and  the  eleven  presentation  points  of  Section  9.1,  we 
compared  (9.2.14)  and  (9.2.15)  against  (9.2.12)  and  obtained  the  statistics  shown 
in  Table  8.  The  statistics  are  given  for  the  results  of  (9.2.12)  and  for  the 
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differences:  [equation  (9.2.14)  or  (9.2.15)]  minus  [result  of  (9.2.12)].  We  should 
note  that  in  view  of  the  results  of  Table  7  we  used  the  integrated  kernel 
evaluation  for  H0  =  3.5  km,  and  simply  used  the  center  point  kernel  evaluation 
for  H0  =  8.5  and  H0  =  28.5  km  to  implement  (9.2.14)  and  (9.2.15). 


Table  8.  Statistics  of  Agtu0  ,  and  of  the  Differences  ((Agt,)^^  -  Agl|}0)  and 
((4gt2)DHo  -  AgfHo),  in  mgals.  See  equations  (9.2.14),  (9.2.15),  and 
(9.2.12). 


Altitude  of 
Upward 
Continuation 

(i) 

A«Ho 

(9.2.12) 

(2) 

(aS»)Bo  "  aIho 

(9.2. 14)-(9.2. 12) 

(3) 

(Ag|)Ho_  a«Ho 
(9.2.15)-(9.2.12) 

H0=28.5  km 

Mean=-3. 11 

-0.63 

-0.02 

S.D.=*4.74 

*0.51 

*0.24 

RMS  =  5.49 

0.80 

0.23 

8.5  km 

-3.71 

-0.64 

0.06 

*11.24 

*0.62 

*0.58 

11.34 

0.87 

0.55 

3.5  km 

-3.39 

-1.58 

*15.57 

*1.22 

*1.04 

15.23 

1.30 

0.99 

The  small  mean  differences  in  Table  8,  column  3  imply  that  numerically: 


M{4gHo>  =  M[(AghDHo}  .  (9.2.16) 


where  M  denotes  the  straight  averaging  operator,  operating  on  the  test  point 
samples.  Therefore,  for  the  difference  (let  us  denote  this  difference  by  Am) 
between  the  third  terms  (9.2.12)  and  (9.2.13)  we  have  the  approximation: 


M{Ain}  =  M{(Agt)j}o  _  4gto}  a  M{(Agt)jj0  -  (AgbBo)  •  (9.2.17) 
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By  virtue  of  (9.2.15b), 


M{(4gt)H2  -  U«2)Ho)  =  M<-(tc  -  t£s)Ho} 


(9.2.18) 


Combining  ((9.2.17)  and  (9.2.18)  we  finally  have  the  approximate  mean  difference 
between  the  third  terms  (ig^^Ho  and  4S^Ho* 


M{Ajij}  =  M{(4gt)Ho  “  4gHo)  *  M{-(tc  -  tcs)n0} 


(9.2.19) 


Combining  (9.2.19)  and  (9.2.8)  we  have: 


M{4IIl}  a  ~M{4l} «  o'" 

M{Aj  +  Ajjj}  a  0.  (9.2.20) 


Prom  the  discussions  related  to  Table  7  we  also  learned  that  M  ( A ji  }  *  0,  and 
therefore  the  mean  of  the  total  difference  A  =  Aj  +  Ajj  +  Ajjj  ia  small,  i.e., 


M{A}  =  M{Aj  +  Ajj  +  Ajjj}  a  0 


(9.2.21) 


Equation  (9.2.21)  repeats  the  results  of  Tables  3,  4,  and  5,  namely,  that  the 
mean  difference  between  the  indirect  method,  on  the  one  hand,  and  the  direct 
method  that  uses  terrain  corrected  surface  anomalies  (Aga  +  tc),  on  the  other 
hand,  is  small. 


9.2.3  Sensitivity  of  Anomaly  Fields  to  Changes  in  Upward  Continuation 
Distance  Hn 

As  mentioned  in  Section  3.3  (A)  there  is  an  uncertainty  as  to  what  value 
of  upward  continuation  distance  H0  to  use  in  the  Poisson  integral.  Recall  that 
H0  is  theoretically  the  vertical  distance  between  the  computation  point  P  in 
space  and  the  level  surface  to  which  the  anomaly  data  are  assumed  to  refer. 
The  uncertainty  in  H0  is  caused  by  the  fact  that  the  given  surface  anomaly 
data  refer  to  a  varying  level  surface,  rather  than  to  a  single  level  surface. 
Even  in  the  case  when  the  terrain  correction  is  applied  to  implement  an 
approximate  data  reduction  to  a  level  surface,  the  position  of  the  final 
reference  level  is  uncertain  (see  Section  3.3  (A)).  In  our  final  procedures  we 
simply  assume  that  the  reference  level  coincides  with  the  mean  elevation 
surface  in  the  area  covered  by  our  anomaly  data.  The  error  in  upward 
continuation  distance  H0  that  results  from  this  assumption  is  expected  to  be 
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related  to  the  deviation  of  the  actual  topography  from  the  mean  elevations 
surface.  In  this  section  we  give  a  feeling  for  the  sensitivity  of  our  results  to 
the  choice  of  value  for  H0. 

The  mean  elevation  in  our  New  Mexico  study  area  is  about  1.5  km.  For 
sensitivity  analysis  we  compared  upward  continuation  results  for  the  case  when 
the  reference  level  for  the  data  is  assumed  to  be  at  the  1.5  km  level,  and  for 
the  case  when  this  level  is  assumed  to  be  at  the  geoid  (i.e.,  the  0  km-levei). 
The  difference  between  the  two  cases  is,  therefore,  the  use  of  upward 
continuation  distances  H0  that  differ  by  1.5  km.  We  performed  our  comparisons 
using  the  eleven  presentation  points  of  the  profile  in  Section  9.1.  The 
statistics  of  the  comparisons,  for  various  values  of  H0  and  various  types  of 
anomaly  fields,  are  shown  in  Table  9. 

Table  9  says,  for  example,  that  in  the  direct  method  that  uses  the  total 
field  (ags  +  tc)DjjOI  (see  (3.3.2)),  an  uncertainty  of  1.5  km  in  H0  for  Hor30  km 
directly  causes  an  uncertainty  of  0.43  mgal  (3.8%  of  computed  value)  in  the 
upward  continued  anomaly.  If  data  reductions  are  used,  as  in  the  indirect 
method,  such  that  only  the  residual  part  (Agr  +  tcs)^H0  is  used  (see  (3.3.14)), 
the  uncertainty  reduces  to  0.19  mgal  (4.2%  of  computed  residual  value).  Table 
9  also  shows  that  the  residual  topographic  field  Ag%0  (see  (3.3.16))  is  the  part 
of  the  total  field  that  is  most  sensitive  to  changes  in  H0,  the  reason  being  that 
agt  contains  the  high  frequency  part  of  the  field.  This  last  sensitivity  does 
not  introduce  any  error  into  the  computations,  since  ig^Ho  has  no  problem 
associated  with  defining  an  "upward  continuation  distance".  Finally,  Table  9 
also  shows  that  the  spherical  harmonic  (long  wavelength)  component  of  the 
field  is  the  least  sensitive  to  altitude  changes. 
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Table  9.  Sensitivity  of  Sample  Profile  Anomalies  to  a  1.5  km  Change  in  Upward 
Continuation  Distance  H0>  for  Various  Values  of  H0  and  Various  Types 
of  Anomaly  Fields.* 


Upward 

Continuation 

Distance 

Total 

Field 

(Aga  +  tc)jjo 

(3.3.2) 

Residual 

Field 

(AgS  +  tc)Ho 

(3.3.14) 

Spherical 

Haraonic 

Field 

A«Ho 

(3.3.15) 

Residual 

Topographic 

Field 

(3.3.16) 

H0  =  30  km 

11.32 

4.49 

9.65 

5.49 

0.43 

0.19 

0.24 

0.27 

3.8% 

4.2% 

2.5% 

4.9% 

H0  =  10  km 

19.92 

8.02 

13.35 

11.34 

1.35 

0.51 

0.37 

1.14 

6.8% 

6.4% 

2.8% 

10.1% 

H0  =  5  km 

23.70 

9.41 

14.59 

15.23 

2.15 

0.73 

0.41 

2.29 

9.1% 

7.8% 

2.8% 

15.0% 
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1.3  Fourier  Computations 

From  the  numerical  point  of  view  it  is  very  attractive  to  use  frequency 
iomain  processing.  The  theoretical  equations  are  summarized  in  Chapter  2.4. 
lere  we  will  show  how  the  Fourier  Transform  principle  can  be  utilized  in 
practice  to  upward  continue  the  gravity  anomaly  field  from  one  level  surface  to 
mother.  It  should  be  stressed  here  that  the  Fourier  technique  requires 
exactly  the  same  assumptions  as  the  flat  Earth  Poisson’s  integral.  After  all 
hese  two  procedures  give  the  unique  solutions  to  the  same  Dirichlet  problem 
'or  half  space  -  one  in  spatial  the  other  in  frequency  domain.  Fourier 
echnique  is  much  cheaper  in  computational  stage  but  requires  the  special  care 
n  controlling  the  edge  effects.  As  our  balloon  experiment  shows,  these  effects 
ire  negligible  far  from  the  edges  where  the  upward  continued  value  was 
:omputed. 

In  practice  the  evaluation  of  upward  continuation  operator  is  done  by 
r-'ans  of  digital  Fourier  transformation  defined  by  (Robinson  and  Silvia,  1981, 
:h.  3.3).  (Notice  that  we  define  the  forward  transform  with  a  +  sign  in  the 
exponent  which  is  a  common  practice  when  working  with  spatial  signals). 


Nx-1  Nv-1 


F(mx,  my)  =  ^  f (nx,  ny)  e 


i2n 


f  Ihd2x  ,  Hy 
l  Nv  N 


i  lylTly 

y 


(9.3. la) 


nx-0  ny=0 


Nx-1  Nv-1 


1  ~ j ~ i 2 7i f x  +HyHly]  (9.3.  ]l,) 

f(nx,  ny)  =  ^  ^  ^  y  F(mx,  my)  e  Nx  Ny  > 

mx=0  my=0 


nx ,  mx  J  { 0 ,  1 ,  .  . .  ,  Nx  1 } 
ny>  my  e{0,  1,  ...,  Ny-1} 


» 


> 


i' 


If  we  choose  Nx,  Ny  to  be  powers  of  2  we  can  use  the  Fast  Fourier  Transform  ^ 

algorithm  to  evaluate  (9.3.1). 

The  integers  nx,  ny,  mx,  my  are  defined  by: 

y  r  ny  ay  , 

ky  =  my  4fy  (9.3.2) 

4fy  =  1 / ( Ny  ay) 

> 


x  -  nx  ax 
kx  mx  a  f x 
afx  =  [/ ( Nx  ax) 


> 
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where  the  spatial  increments  are  related  to  angular  increments  by 


Ax  -  R  COS  ♦  AX 
Ay  =  R  A$ 


(E-W  direction) 
(N-S  direction) 


(9.3.3) 


Here  the  regular  spatial  grid  intervals  Ax  and  Ay  are  computed  using  the  mean 
latitude  ♦=?=32.5*  of  the  data  location.  The  original  data  are  considered  here 
as  the  discrete  point-measurements  of  the  gravity  anomaly  field  of  the  surface 
of  the  earth  using  uniform  angular  spacing  A<>rAXr5’.  This  uniform  spacing  on 
the  sphere  under  transformation  (9.3.3)  becomes  strictly  speaking  non-uniform 
on  the  plane  due  to  meridian  convergence.  Using  a  fixed  value  of  ♦  in  (9.3.3) 
introduces  some  distortions  to  the  original  spacial  distributions  of  the  data. 
Our  study  indicates  that  this  distortion  produces  an  essentially  negligible 
effect  on  the  results  due  to  the  use  of  rather  fine  5’x5’  grid  defined  on  a 
relatively  small  portion  of  the  sphere  (7*x9°  block).  In  other  words  the  flat 
earth  approximation  being  valid  for  the  original  Poisson  space  processing  •' 
data  is  also  valid  for  its  frequency  implementation,  at  least  for  small  regions  as 
discussed  in  this  report. 

Summary  of  frequency  domain  digital  upward  continuation  procedures: 

1.  compute  Ax  and  Ay;  (9.3.3) 

2.  transform  the  data  to  the  frequency  domain  as  indicated  by  (2.4.5)  using 
(9.3. la); 

3.  multiply  the  transformed  data  in  the  frequency  domain  by  the  transfer 
function  A  (2.4.6) 

4.  invert  the  result  back  to  space  domain  as  indicated  by  (2.4.7)  using 
(9.3.1b). 


The  above  procedure  was  applied  to  our  7*x9*  New  Mexico  study  area,  to 
upward  continue  the  entire  84x109  grid  of  5’x5’  mean  Faye  anomalies  (see 
Section  8.7).  The  uplifted  grid  for  an  upward  continuation  distance  of  28.5  km 
is  contoured  in  Figure  16,  with  a  contour  interval  of  5  mgals.  This  figure 
should  be  compared  with  Figure  17  in  which  the  corresponding  grid  from  the 
Rapp-180  field  is  contoured.  The  long  wave  agreement  between  Figures  16  and 
17  is  evident,  with  Figure  16  expectedly  showing  more  detail  because  of  the 
use  of  high  frequency  terrestrial  data. 

For  the  Fourier  upward  continuation  we  used  only  a  simplified  procedure 
including  no  windowing  of  the  data  or  any  other  regularization  routine  in 
frequency  domain  (i.e.,  no  specialized  filters  used).  We  found  the  results 
produced  by  this  simplified  procedure  in  satisfactory  agreement  with  the  space 
domain  processing  of  original  Ag+tc  field  by  means  of  Poisson’s  integral.  The 
difference  in  results  on  the  order  of  0.15  mgal  at  30  km  flight  altitude  can  be 
associated  with  the  differences  iri  numerical  implementation.  We  employed  this 
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‘iguro  16.  Upward  Continued  Faye  Anomalies  Using  Fast  Fourier  Upward 
Continuation,  New  Mexico  Balloon  Gravity  Test  Site. 

C.I.  =  5  mgals;  Upward  Continuation  Distance  =  28.5  km. 


LRT I TUDE 


'dure  using  the  Fast  Fourier  Transform  routine  m  Klifb'd  f r< >:n  > 'laerl.o,.! 

,  p.  12).  The  technique  was  tested  using  the  original  terrain  ■■■urn-. -fed 
air  anomaly  field.  The  results  were:  compared  along  K-W  test  profile  with 
•esults  of  the  Poisson’s  integral  operating  on  the  same  field. 

The  RMS  difference  over  the  test  profile  between  those  two  .results  wn 
mgnl  for  the  upward  continuation  distance  of  28. f>  krn.  For  the  upward 
nuation  distance  of  8.5  km  the  RMS  difference  along  the  same  test  profile 
3.32  tngal. 

The  magnitude  of  the  above  differences  is  in  fact  on  the  order  of 
rences  between  the  indirect  and  direct  space  domain  procedures  already 
issed.  This  indicates  that  the  Fourier  method  can  prove  to  be  competitive 
the  one-step  space  domain  processing  (direct  method)  with  respect,  to 
-ncy  of  the  results.  It  can  be  recommended  for  fast  processing  of  large 
nes  of  gridded  data  to  produce  the  image  of  the  field  (in  a  gridded  form) 
ly  level  surface  above  the  data  surface.  To  get  the  values  at.  any  location 
de  the  grid  point  the  interpolation  has  to  be  performed  which  is  a  valid 
ne  provided  the  original  data  were  gridded  in  accordance  with  the 
ling  theorem  (Robinson  and  M.T.  Silvia,  1981,  eh.  2.6). 

)n  the  theoretical  side  it  provides  an  elegant  uniform  treatment  of  all 
jencies  present  in  the  original  signal.  On  the  practical  side,  it  provides  a 
numerical  procedure  to  transform  the  data  upward  from  one  level  surface 
nother.  As  far  as  accuracy  is  concerned  the  edge  effects  (inherent  in 
ier  techniques)  should  be  treated  with  care  by  appropriate  windowing 
ne  or  (if  possible)  the  region  with  data  should  cover  enough  grounds 
Hireling  the  point  of  interest. 


^plications  to  Balloon  Gravity  Project 

to  theory  outlined  in  this  report  was  applied  to  the  Balloon  Gravity 
:t  coordinated  by  the  Air  Force  Geophysics  Laboratory,  Bedford, 
husotts.  The  experiment  took  place  in  New  Mexico  and  was  designed  to 
he  theory,  procedures  and  instruments  used  for  both  the  measurements 
it1  prediction  of  gravity  in  space. 

he  comparisons  between  the  observed  and  predicted  gravity  will  give  in 
into  the  accuracy  and  performance  of  the  theories  and  techniques  of 
y  recovery  in  space  which  are  in  the  operational  stage  today.  It  is  also 
.  of  the  accuracy  and  performance  of  the  balloon-borne  instruments  and 
iques  that  are  used  today  for  the  measurements  of  the  external  gravity 

1  this  section  some  technicalities  of  the  actual  gravity  prediction 
dure  used  for  the  Balloon  Project  are  given.  The  method  we  chose  for 
nal  application  is  the  indirect  method  described  in  Section  3.3.  To  repeat, 
y,  the  rationale  behind  this  method  is  to  extract  the  high  frequency 
ion  in  the  original  surface  gravity.  Here  by  high  frequency  variation  we 
■stand  the  topographic  effects  (irregular  geometry  of  the  terrain  plus  tc 
dion).  As  we  have  seen  in  the  previous  chapters  this  high  frequency 
merit  of  the  gravity  signal  cannot  be  conceptually  nor  practically 
rted  downward  to  the  Poisson’s  spherical  geometry  so  we  did  not  try  to 
lis.  Instead  we  immediately  upward  continued  this  high  frequency 
merit  (right  from  the  original  surface  where  it  was  defined)  up  to  the 
m's  attitude  using  partially  the  ideas  of  equivalent  source  method  as 
minted  in  the  prism  integration  of  topographic  effects. 

f ter  the  removal  of  this  shor  t-wavelength  variation,  the  remaining  regular 
>n  of  the  gravity  field  was  converted  downward  (formally)  to  the  Poisson’s 
■ical  geometry  (see  Section  3.3).  Then  the  long-wavelength  or  global 
m  of  this  signal  is  upward  continued  in  the  frequency  domain  (using 
•ical  harmonics)  and  the  mid-frequency  portion  is  upward  continued  by 
lg  the  usual  Diriehlet.  problem  for  half  space  (in  planar  approximation)  by 
<  f  Poisson’s  integral. 

'  !■•■  a  ■'  ual  ■  at  put  at  ion  we  used  the  FORTH  V-.'  imelomeut.V  ion  of  the 
m’s  integral  procedure  as  described  in  (Rapp,  February  19B6).  The 
to  was  run  using  the  gridded  mid-frequency  portion  of  the  gravity 
»l>  signal  as  discussed  in  Section  3.3. 


I’reparat  ion  of  the  Balloon  Tracking  Data  for  Upward  Continuation 
[  ’p  e,  >rj  |j  r«' 

he  balloon  tracking  data  were  sent  on  the  magnetic  tape  dated  4  May 
Hie  tape  name  is  DUCKY  fa  Radar  Tracking  Data.  The  tape'  contained 
files,  with  the  float  portion  of  the  experiment  contained  as  part  of 
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The  positions  of  the  balloon  during  the  experiment  were  provided  in  the 
form  of  geodetic  coordinates  (♦,  X,  h)  with  respect  to  the  WGS72  ellipsoid  (see 
Table  10). 

The  original  tape  contained  462,200  positions.  The  original  records  were 
parameterized  by  the  time  of  measurements.  On  Figure  18  we  show  the  spatial 
location  of  the  balloon’s  trajectory.  The  average  spacing  of  original  tracking 
data  for  the  float  portion  of  the  experiment  was  approximately  0.6  m  on  the 
ground.  At  the  first  step  we  reduced  the  number  of  data-points  spatially  by 
choosing  only  the  clusters  of  10  original  records  spaced  every  30"  angular 
distance  apart  from  each  other.  During  this  step  the  false  data  records  were 
rejected  and  the  remaining  set  containing  3160  records  was  checked  for 
blunders.  The  false  data  records  on  the  original  tape  occurred  at  the  end  of 
each  file  due  to  technical  reasons. 

Since  we  were  primarily  interested  in  the  flight  altitude  portion  of  the 
experiment  (29-30  km  altitude)  it  was  sufficient  to  select  only  1  data  point 
every  2’  spacing  (in  angular  distance)  along  the  track.  The  2'  spacing  is 
sufficient  for  all  interpolation  purposes  at  flight  elevation  because  the 
anomalous  gravity  field  is  already  very  smooth  at  that  altitude  (see  Figure  15a, 
for  example).  The  resultant  data-set  contained  31  values  equally  spaced  in  2’ 
intervals  covering  the  flight  portion  of  the  balloon’s  trajectory,  that  is  the 
portion  of  the  original  data  tape  (file  #3)  which  falls  in  the  time  interval 
<57445.95,  66955.95>.  The  time  is  UTC  time  in  seconds.  In  the  computational 
stage  of  this  project  we  used  the  time  only  as  a  convenient  parameter  to  locate 
the  data  on  tape. 


10.2  Reference  Systems  Used  in  the  Balloon  Project 

In  this  section  we  give  the  summary  of  reference  systems  and  conversions 
used  at  the  stage  of  data  preparation  for  the  balloon  project.  The  geodetic 
coordinates  (♦,  X,  h)  of  the  balloon  positions  provided  on  balloon  tape  were 
given  with  respect  to  WGS72.  The  terrestrial  gravity  anomaly  field  that  we 
used  in  this  project  was  given  in  the  GRS67  system.  The  spherical  harmonic 
expansion  of  the  global  gravity  field  up  to  degree  180  was  assumed  to  refer  to 
GRS80  (Rapp,  1981)  (referred  to  as  Rapp  180  field).  The  ellipsoid  parameters 
for  the  reference  systems  used  are  given  in  Table  10. 


Ellipsoid 

a  [m] 

1/f 

WGS72 

6378135 

298.26 

GRS67 

6378160 

298.2471674273 

GRS80 

6378137 

298.257222101 

Table  10.  Parameters  of  the  reference  ellipsoids  involved  in  the  datu 
preparation  for  the  balloon  project. 


Now  we  will  summarize  the  conversions  used  for  upward  continuation  of 
gravity  anomalies.  Since  the  original  gravity  anomaly  field  that  we  used  in  our 
project  refers  to  GRS67  ellipsoid  we  decided  to  convert  the  given  geometric 
heights  hvvGS72  that  refer  to  the  WGS72  ellipsoid  to  the  normal  heights  H*qRS67 
that  refer  to  the  GRS67  ellipsoid. 


geop 


spherop 


geoid 


ellipsoid 


Figure  19.  The  Spatial  Relationship  Between  Selected  Gravimetric  Quantities 


From  Figure  19  we  have  the  relation: 


H*  =  ii  +  (n  -  ()  _  (10.2.1) 

GRS67  V  WGRS67 

where  (  is  the  geop-spherop  separation  or  height  anomaly  at  bnlloon’s 
altitude. 

Instead  of  the  difference  (N-{)grs67  we  will  actually  use  (N-OGRS80  in  the 
approximate  relation: 


H*GRS67  3  H  +  (N  f)GRS80 


(10.2.2) 
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For  the  orthometric  height  II  appearing  in  (10.2.2)  we  can  write  the  relation 
(see  Figure  19): 


H  "  hGRS80  NGRS80 


(10.2.3) 


To  get  h(-jRS80  from  a  given  hwGS72  we  implement  the  well  known  formulas  of 
geometric  geodesy  (Rapp,  1984,  Geometric  Geodesy  Notes,  Vol.  1,  pp.  121-122)  in 
the  following  procedure 


input: 

compute: 

convert : 


compute: 


^WGS72 ’  * 

ZWGS72  =  W1  -  e2)  +  h)  sin  ♦  using  WGS72 
constants 

ZGRS80  =  ZWGS72  +  42  ZWGS72 

here  we  set  the  origin  shift  AZ  to  zero  so  that 
the  ellipsoid  WGS72  and  GRS80  have  a  common  center 

h /-.rioor,  "  -  N  +  e2N  using  GRS80  constants 

GRS80  sin  ♦ 


(10.2.4) 


where  N  =  n//l-e2sin24>  ,  e2  =  2f  -  f2. 


In  practice,  instead  of  converting  each  data  point  separately  using  the 
above  formulas,  we  applied  a  single  common  constant  Ah  =  -1.66  m  to  each 
height  given  in  the  WGS72  system: 


hGRS80  hWCS72  ‘  4h 


(10.2.5) 


where  the  particular  value  of  Ah  -  hQRS80-hWGS72  =  -1.66  m  has  been 

evaluated  by  the  sequence  of  equations  (10.2.4)  (starting  from  a  nominal 
hwos72  and  computing  the  equivalent  hcRS80)  using  the  mean  latitude  of  the 
balloon’s  trajectory  (32.40*)  and  the  mean  altitude  of  the  flight  section  of  the 
experiment  (30  km). 

Using  (10.2.3)  and  (10.2.5)  we  can  write  (10.2.2)  in  the  final  form: 


GRS67 


s  ( 


"WGS72 


Ah) 


'  GRS80 


(10.2.6) 


In  the  actual  implementation  of  equation  (10.2.6)  we  used  the  Rapp  180  field 
(based  on  the  spherical  harmonic  expansion  of  gravity  potential  up  to  degree 
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and  order  180)  to  generate  the  required  values  of  height  anomalies  {‘GRS80- 
Due  to  the  smoothness  of  the  height  anomaly  field  at  flight  altitude  a  single 
constant  mean  value  of  {GRS80  =  -23  m  (*0.26  m)  was  used  in  formula  (10.2.6). 

The  normal  heights  (10.2.6)  were  subsequently  used  in  the  FORTRAN 
subroutine  for  Normal  Field  Computations  (Rapp,  1966)  (see  also  Section  2.6)  to 
generate  normal  gravity  at  the  exact  balloon  positions.  The  routine  was  run 
using  GRS67  constants. 


10.3  Upward  Continuation  Results 

The  upward  continuation  procedure  of  gravity  anomaly  described  in  this 
report  was  actually  applied  only  to  the  31  points  selected  in  Section  10.1.  The 
results  are  shown  in  Table  11.  The  position  of  each  point  is  defined  in 
columns  1  to  4:  (1)  time  tag  in  the  tracking  record;  (2)  latitude;  (3)  longitude; 

and  (4)  height  above  the  WGS72  ellipsoid.  Columns  5  to  7  give  the  upward 
continued  anomaly  contribution  from:  (5)  residual  defined  Bouguer  anomalies; 
(6)  spherical  harmonic  anomalies;  and  (7)  topographic  anomalies.  Columns  5,  6, 
7  were  summed  to  form  column  8,  which  is  the  total  upward  continued  anomaly 
from  the  indirect  method  of  upward  continuation.  Column  9  gives  the  normal 
gravity  to  be  added  to  the  upward  continued  anomaly  to  produce  the  first 
model  for  measured  gravity  at  the  space  point.  Columns  (10)  and  (11)  give  the 
upward  continued  anomalies  resulting  from  the  direct  method  of  upward 
continuation  of  surface  anomalies  using  terrain-uncorrected  and  terrain- 
corrected  surface  anomalies,  respectively.  Columns  (12)  and  (13)  give  the 
errors  of  columns  (10)  and  (11)  relative  to  the  expectedly  more  rigorous 
indirect  method  (column  8)  of  upward  continuation.  The  mean  error  and 
standard  deviation  of  errors  are  as  follows: 


Error  of  Direct  Method 
without  use  of  tc  (col.  12) 


Error  of  Direct  Method 
with  use  of  tc  (col.  13) 


mean  error 


-0.51  mgal 


0. 14  mgal 


s.d.  error 


*0.35  mgal 


*0. 18  mgal 


In  agreement  with  earlier  results,  we  see  an  improvement  with  the  use  of 
terrain  correction  in  the  direct  method. 


10.4  Interpolation  of  the  Results  at  Altitude 

In  the  last  step  the  upward  continued  values  were  interpolated  at  all 
original  data-records  on  the  balloon  tape  that  fell  in  the  time  interval  covering 
the  flight  portion  of  the  experiment  (190,201  data-points).  Only  results  from 
the  indirect  upward  continuation  method  were  used  (Table  11,  Col.  8.).  The 
actual  interpolation  was  done  using  the  one-dimensional  cubic  spline  routine 
(subroutines  SPLINE  and  SEVAL,  Forsythe  et  al.,  1977)  applied  to  the  flight 
sequence  which  has  been  parameterized  with  longitude  only. 
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Tublc  11.  Gravity  Anomaly  and  Normal  Gravity  Upward  Continuation  Results  at  31  Balloon  Points,  New  Mexico 
Area.  Units:  mgals. 
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To  check  the  internal  accuracy  of  this  interpolation  process  we  removed 
every  other  record  out  from  the  31  "master"  points  and  treated  the  remaining 
16  points  as  the  new  knots  in  the  spline  interpolation  procedure.  In  this  way 
we  interpolated  values  at  the  removed  knots  and  compared  with  the  "true" 
values  obtained  directly  from  the  upward  continuation  process.  The 
differences  were  on  the  order  of  0.01  mgal  in  gravity  anomaly. 

For  the  normal  gravity  computation  the  chosen  31  grid  points  were  not 
dense  enough  to  properly  recover  the  normal  gravity  at  flight  altitude,  mostly 
due  to  the  vertical  gradient  of  normal  gravity.  It  was  evident  especially  at 
both  ends  of  the  flight  portion  of  the  experiment  where  the  balloon  changes 
altitude  very  rapidly.  Therefore  we  decided  not  to  interpolate  normal  gravity 
but  to  compute  it  rigorously  at  every  observational  point  (see  Section  2.6). 


10.5  Propagation  of  Positional  Errors 

On  the  original  balloon  tape  the  positional  accuracy  of  the  balloon  in  all 
three  directions  were  provided  with  each  data-record.  The  average  accuracy 
of  the  31  points  selected  for  actual  upward  continuation  was  on  the  order  of 
2  m  in  all  x,  y  and  z  coordinates. 

Considering  that  the  largest  horizontal  gradient  of  the  actual  computed 
anomalous  gravity  field  along  the  flight  portion  of  the  balloon  trajectory  is 
about  0.0006  mgal/m  (this  is  the  actual  gradient  of  the  computed  results  for 
the  balloon  project),  we  assess  the  maximal  error  due  to  the  2  m  uncertainty 
in  the  horizontal  position  of  the  balloon  to  be  0.0012  mgal  in  gravity  anomaly 
and  0.0011  mgal  in  normal  gravity  (for  normal  gravity  computation,  only  the 
uncertainty  in  N-S  direction  must  be  considered). 

The  uncertainties  in  altitudes  will  also  produce  uncertainties  in  the 
computed  components  of  the  gravity  field.  From  Table  9  we  learn  that  the 
vertical  gradient  of  the  actual  Ag  field  at  30  km  flight  altitude  is  about  0.0003 
mgal/m.  Therefore  the  2  m  error  in  the  balloon  vertical  position  will  show  up 
as  the  uncertainty  of  about  0.0006  mgal  in  the  computed  Ag  field.  Similarly 
the  vertical  gradient  of  normal  gravity  of  0.3086  mgal/m  (Rapp  1982,  p.  8)  will 
cause  the  uncertainty  in  the  computed  normal  gravity  of  about  0.6  mgal  (due  to 
the  2  m  uncertainty  in  the  altitude). 

Also,  the  error  in  geoid  undulation  (more  precisely  in  height  anomalies) 
used  in  the  data  reduction  propagates  directly  to  the  uncertainties  in  the 
clearance  of  the  balloon  above  the  datum  surface.  In  our  project  we  used  the 
height  anomalies  field  implied  by  the  set  of  potential  coefficients  up  to  degree 
180  computed  by  (Rapp,  Dec.  1981)  and  referred  to  as  the  Rapp  180  field.  The 
height  anomalies  were  needed  in  data  reduction  process  to  convert  the  WGS72 
geometric  heights  given  as  the  data  to  the  normal  heights  in  GRS67  (See 
Section  10.2). 

The  Rapp- 180  field  used  in  the  data  reduction  process  gives  the 
undulations  (height  anomalies)  with  the  accuracy  on  the  order  of  *1  m  (Rapp, 
Dec.  1981,  p.  31).  For  our  actual  sub-balloon  trace  in  New  Mexico  we  decided  it 
is  sufficient  to  represent  the  Rapp-180  height  anomaly  field  in  this  area  by  a 
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constant  value  of  -23  m  (*0.26  m)  which  is  the  mean  height  anomaly  computed 
for  the  flight  portion  of  the  experiment  only.  We  conclude  that  the  usage  of 
the  Rapp-180  undulation  field  introduces  about  1  m  uncertainty  into  the  balloon 
vertical  position.  Using  the  same  vertical  gradients  stated  above  we  find  the 
0.0003  mgal  error  in  the  computed  Ag  value  and  the  0.3  mgal  error  in  the 
computed  normal  gravity. 

In  conclusion,  we  notice  that  the  positional  errors  in  the  balloon 
coordinates  and  the  height  anomaly  error  do  not  noticeably  affect  the  gravity 
anomaly  computation,  but  they  have  an  effect  on  the  computed  normal  gravity 
and  on  the  predicted  observed  gravity  (mainly  due  to  the  altitude  error). 
Therefore,  the  predicted  observed  gravity  may  be  contaminated  by  at  least  0.6 
mgal  uncertainty  due  to  the  vertical  positional  errors. 


10.6  Other  Sources  of  Errors 


In  Section  6  we  considered  theoretically  how  the  errors  present  in  gravity 
data  propagate  through  the  upward  continuation  integration.  Here  we  can  try 
to  use  some  of  the  results  of  Section  6  to  give  the  rough  evaluation  of  errors 
that  affect  the  actual  procedure. 

The  main  source  of  the  gravity  data  error  probably  comes  from  the 
gridding  procedure  of  the  gravity  material  by  means  of  the  collocation  pre¬ 
diction  of  the  mean  5’x5’  values.  If  we  assume  that  the  mean  values  computed 
by  the  collocation  prediction  procedure  are  contaminated  by  the  error  function 
having  error  variance  of  25  mgal2  (this  number  is  a  rough  estimate  from  the 
formal  errors  output  by  the  collocation  prediction  in  the  area  of  the  balloon 
flight;  see  Figure  12.  central  portion)  and  the  correlation  length  of  about  10 
km  (this  correlation  length  is  unavoidably  due  to  the  gaps  in  the  original 
gravity  data  of  about  this  size  which  cause  the  adjacent  5’x5’  blocks  to  be 
correlated),  then  we  can  use  Table  2  together  with  (eq.  6.2.9)  to  conclude  that 
at  30  km  flight  altitude  the  propagated  effect  is  about  0.9  mgal  (in  standard 
error)  with  55  km  correlation  length.  This  puts  the  limit  of  accuracy  on  our 
entire  procedure.  This  limit  is  due  to  the  quality  of  the  original  gravity  data 
mainly  in  spatial  distribution  and  cannot  be  overcome  by  refinements  in 
procedure  unless  the  geometry  and  quality  of  original  data  are  improved. 

It  should  be  noted  here  that  errors  present  in  the  mean  anomalies  (for 
which  25  mgal2  variance  and  10  km  correlation  length  at  ground  level  was  our 
rough  estimate)  are  due  to  both:  errors  in  the  original  gravity  point-values 
(on  input  to  collocation  prediction  routine)  and  errors  of  interpolation.  The 
error  of  interpolation  comes  from  the  difference  between  the  true  gravity 
anomaly  and  the  gravity  anomaly  model  implemented  by  the  subtraction  of  the 
reference  field  in  order  to  center  the  original  data.  Of  course  where  the 
original  data  are  very  dense  the  computed  mean  values  are  (almost) 
uncorrelated  and  are  affected  mainly  by  the  point-gravity  data  error.  If  this 
error  in  original  point-values  could  be  modelled  by  a  weakly  correlated  noise 
having  correlation  length  shorter  than  5’x5’  blocks  (used  in  prediction)  then 
this  type  of  data  error  would  tend  to  cancel  during  the  computation  of  5’x5’ 
means,  producing  essentially  negligible  effect  at  30  km  altitude.  If  the 
point-data  field  is  sparse  in  some  areas  then  the  predicted  mean  values  will  be 


ffected  by  both  data  error  and  the  interpolation  error.  In  sparse  areas  the 
omputed  mean  values  tend  to  be  correlated  with  each  other  and  so  will  be  the 
interpolation  error.  The  correlation  length  of  this  error  is  dependent  on  the 
lize  of  the  gaps  in  the  original  data  -  the  larger  the  gaps  the  wider  correlated 
trrors  (from  this  type  of  analysis  comes  our  rough  estimate  of  10  km 

;orrelation  length  errors  present  on  output  of  collocation  prediction 

jrocedure).  As  we  learn  in  Section  6  such  widely  correlated  errors  do  not 
ittenuate  very  fast  with  altitude.  For  example,  according  to  our  rough 
tstimate  (see  Table  2)  at  30  km  flight  altitude  the  correlated  errors  produce 

,he  effect  that  can  be  descibed  as  the  distortion  function  which  is  0.9  mgal  in 

amplitude  and  55  km  wide  in  correlation  length.  This  type  of  effect  can  very 
>asily  be  misinterpreted  as  some  sort  of  systematic  error  and  wrongly 
issociated  with  errors  in  the  modeling  of  the  upward  continuation  distance  (see 
Section  9.2.3)  or  the  data  reduction  error.  We  consider  this  type  of  error  as 
Lhe  main  limitation  of  accuracy  of  final  results. 

Another  effect  coi.sidered  in  Section  6  is  the  representation  error  due  to 
Lhe  conceptual  replacement  of  the  true  gravity  field  by  the  step  function 
composed  of  the  flat  patches  over  5’x5’  blocks.  At  30  km  flight  altitude  the 
rough  estimate  of  this  effect  in  gravity  anomaly  is  only  about  0.01  mgal 
(standard  error). 

If  we  assume  tha)  the  effects  considered  here  act  independently  of  each 
other  we  can  sum  th  !  error  variances  of  each  component  to  get  the  rough 
estimate  of  the  total  ?rror  to  be  0.92  mgal  (total  standard  error)  in  gravity 
anomaly  and  about  0.7  mgal  error  in  normal  gravity. 

In  Table  12  we  give  the  summary  of  errors  that  affect  the  actual  results  of 
upward  continued  gravity  anomalies  and  normal  gravity  for  the  Balloon  Project. 
Notice  that  the  indirect  method  is  to  some  extent  free  of  the  truncation  error 
(Section  4)  since  we  carry  up  the  complete  global  long-wavelength  information 
represented  by  the  spherical  harmonic  expansion  of  the  gravity  field  up  to 
degree  180. 

In  order  to  get  a  feeling  of  the  actual  magnitude  of  propagated  errors  it  is 
possible  to  perform  s  me  simplified  numerical  test.  (Rapp,  1966)  suggested  a 
simple  numerical  chec.t  (not  estimate)  on  the  relative  magnitude  of  propagated 
uncorrelated  errors  n  the  upward  continuation  process.  The  trick  is  to 
upward  continue  the  estimates  of  accuracies  of  gravity  material  using  the  same 
computer  program  which  was  used  to  process  the  actual  gravity  data. 

In  a  computer  implementation  the  upward  continuation  integral  (6.1.1)  takes 
the  form  of  summation: 
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(10.6.1) 


This  is  a  weighted  average  operator.  Now,  if  we  know  the  accuracy  m2ag  of 
each  component  Agjj  used  in  the  summation  (10.6.1)  we  can  neglect  the 


jrrelations  between  adjoining  blocks  (i,  j)  and  simply  sum  the  individual  error 
ariances  m*Ag  over  all  blocks  used  according  to  the  formula: 


■  <2^  LL  "’*«  W  (10-6-2) 


'able  12.  Error  Budget  of  the  Gravity  Upward  continuation  for  the  Balloon 
Project 


Effect  at  30  km  Altitude 

Source 

Gravity  Anomaly 

Normal  Gravity 

The  original  gravity  material 
error  and  the  error  in  predicted 
5’x5’  mean  anomalies  (interpola¬ 
tion  error)  (see  Section  6) 
assumed:  standard  error  5  mgals 

correlation  length  10  km 

0.9  mgal  standard 
error 

55  km  correlation 
length 

unaffected 

Errors  in  the  modelling  of  the 
upward  continuation  distance 
(on  the  transition  from  the  true 
earth  to  the  Poisson’s  spherical 
geometry).  (See  Table  9,  col.  3). 
assumed  uncertainty:  1.5  km 

0.2  mgal  uncer¬ 
tainty 

unaffected 

Representation  error 
(See  Section  6) 

0.2  mgal  standard 
error 

unaffected 

Errors  in  the  balloon  position 
(See  Section  10.5) 

horizontal  errors  of  2  m 

0.0012  mgal 
maximal  error 

0.0011  mgal 
uncertainty 

vertical  errors  of  2  m 

0.0006  mgal 
uncertainty 

0.6  mgal  un¬ 
certainty 

Height  anomaly  error 
(See  Section  10.5) 

Estimated  uncertainty  of  1  m 

0.003  mgal 
uncertainty 

0.3  mgal 
uncertainty 
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is  the  quadrature  formula  similar  to  (10.6.1)  but  with  the  squared  kernel, 
•efore,  we  used  the  actual  upward  continuation  program  described  in  (Rapp 
)  to  upward  continue  the  uncorrelated  errors  in  the  gravity  material  due 
Lhe  prediction  of  aean  values  on  the  5'x5’  grid,  using  (10.6.2).  For 
{ramming  details  see  (Rapp,  1966).  Computation  was  performed  exactly  at 
balloon  locations  giving  the  accuracies  on  the  order  of  0.2  mgal  for  the 
it  portion  of  the  balloon  trajectory  (30  km  upward  continuation  distance), 
xi  this  single  experiment  we  observe  the  attenuation  of  errors  from  about  4 
Is  (on  average)  at  the  ground  level  to  about  0.2  ragals  at  30  km  flight 
.ude. 

It  is  important  to  realize  that  this  numerical  experiment  is  valid  only  for 
specific  type  of  error  (namely  errors  due  to  prediction  of  the  mean  gravity 
maly  values  on  the  grid,  see  Section  8.3  for  details)  assuming  the  error  is 
orrelated.  Also  the  arbitrary  scaling  of  error  variances  on  the  ground 
ild  produce  the  respective  rescaling  of  upward  continued  variances  of  flight 
ration.  In  that  sense  only  the  relative  degree  of  attenuation  is  a  meaningful 
come  of  that  numerical  check. 


Summary  and  Conclusion 

We  have  presented  operational  procedures  for  the  upward  continuation  of 
ity  anomalies  given  on  the  surface  of  the  earth.  The  main  conceptual 
culty  is  that  operationally  available  anomalies  are  referred  to  the  earth’s 
ace,  which  is  not  an  equip’otential  surface.  These  surface  anomalies  can  be 
ard  continued  using  discrete  estimation  procedures  such  as  collocation  or  a 
hammar-type  of  approach,  but  in  this  report  we  have  avoided  such 
niques  because  of  the  expensive  matrix  inversions  that  they  require. 

Instead,  as  stated  in  Section  1  we  used  collocation  or  ly  in  a  preliminary 
i,  to  predict  an  optimal  set  of  mean  anomaly  values  from  the  available 
gularly  distributed  point  anomaly  data.  This  application  of  collocation  is 
rationally  feasible  because  in  contrast  to  the  prediction  of  upward 
Linued  values,  the  prediction  of  mean  anomalies  requires  information  to  be 
jrted  only  from  a  small  number  of  data  points  around  the  computation  block, 
sr  obtaining  a  complete  set  of  mean  anomalies  over  rectangular  blocks  we 
i  to  a  continuous  upward  continuation  problem,  in  which  it  is  assumed  that 
;very  point  on  the  earth’s  surface  we  know  the  gravity  anomaly  function,  as 
resented  by  the  mean  values. 

The  upward  continuation  of  a  continuous  gravity  anomaly  function  given 
the  (non-level)  surface  of  the  earth  is  by  no  means  a  simple  problem.  The 
plest  conceptualization  of  a  solution  is  by  means  of  Taylor  series  expansion, 
which  the  surface  anomalies  are  first  used  to  derive  anomalies  on  a  level 
face  using  the  vertical  gradients  of  the  anomaly  field.  Once  the  level 
face  anomalies  are  known  classical  Poisson  integration  yields  a  solution  to 
upward  continuation  problem  with  relative  accuracy  on  the  order  of  the 
th’s  flattening.  However,  for  rough  anomaly  fields  the  computation  of 
•tical  gradients  required  for  data  reduction  to  a  level  surface,  itself  requires 
:h  density  and  accuracy  of  data  that  is  not  usually  available  in  practice  (see 
>  (1980),  for  example).  Moreover,  for  such  fields  downward  continuation  of 
•face  values  cannot  be  expected  to  be  regular.  Therefore,  as  reasoned  out 
Section  3  we  have  resorted  to  the  so-called  indirect  method  of  upward 
itinuation. 

The  most  important  feature  of  the  indirect  method  is  the  extraction  and 
>arate  modeling  of  the  high  frequency  irregularities  of  the  original  gravity 
>mal y  signal.  This  high  frequency  component  is  due  to  shal  ow  topographic 
gses,  and  can  be  modeled  by  directly  integrating  the  gravitational  effects  of 
pse  masses  without  need  for  any  sort  of  data  reduction  to  a  level  surface, 
srationally,  the  topographic  effects  on  gravity  anomalies  at  altitude  can  be 
nputed  by  prism  integration  as  stated  in  Section  2.5  with  an  operational 
)gram  mentioned  in  Section  7,  while  on  the  earth's  surface  the  topographic 
ect  on  gravity  anomaly  at  a  point  is  conveniently  formed  as  the  sum  of  a 
puguer  plate"  effect  and  a  "terrain  correction"  effect  (tee  (3.3.11)).  In  this 
port  we  modeled  the  high  frequency  component  of  the  f -dimensional  gravity 
jmaly  field  using  as  "equivalent  sources”  the  positive  and  negative 
pographic  masses  of  assumed  density  2.67g/cm3,  lying  between  the  actual 
•ography  and  a  reference  topography  to  spherical  harmonic  degree  and 
ler  180.  With  the  extraction  and  separate  modeling  of  the  high  frequency 
jmaly  signal  we  circumvent  the  major  difficulties  associated  with  the 
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reduction  of  surface  data  to  a  level  surface,  these  difficulties  being  precis  ty 
due  to  the  high  frequency  irregularities  of  the  field. 


The  residual  field  left  after  removal  of  topographic  effects  is  much 
smoother  than  the  original  field.  Prom  this  residual  field  we  decided  to 
further  remove  and  separately  model  the  low  frequency  component  using  the 
Rapp-180  (1981)  field.  This  was  done  in  order  to  formally  free  the  upward 
continuation  from  truncation  error  caused  by  neglect  of  remote  zone  data.  In 
order  to  remove  the  effect  of  the  Rapp-180  field  from  the  surface  data,  the 
data  points  were  basically  taken  at  their  actual  horizontal  positions,  but  an 
assumption  had  to  be  made  that  the  data  points  all  lie  on  a  common  level 
surface  and  not  in  their  actual  vertical  positions.  The  assumption  was 
necessary  to  keep  evaluation  time  for  the  Rapp- 180  field  reasonable.  The 
assumption  seems  justified  because  the  vertical  gradient  of  a  180-field  is 
expected  to  be  small,  but  this  point  can  be  further  studied  (see,  for  example 
Table  9).  The  Rapp-180  field  was  evaluated  at  ground  level  using  a  program 
for  fast  generation  on  a  grid,  while  at  isolated  computation  points  at  altitude 
another  program  suited  for  single  point  computations  was  used  (see  Section  7). 

The  medium  frequency  residual  field,  left  after  removing  both  the  high 
frequency  topographic  effects  and  the  low  frequency  spherical  harmonic  field, 
was  then  modeled  by  the  Poisson  integral.  Since  the  data  points  were  still 
located  on  the  earth’s  surface,  a  data  reduction  to  a  level  surface  was  still 
called  for.  However,  since  the  residual  field  is  much  smoother  than  the 
original  field,  an  approximate  reduction  can  be  used.  To  do  this  a 
long-wavelength  form  of  the  terrain  correction,  namely,  the  terrain  correction 
tcs  of  an  expansion  of  the  topography  to  cegree  180,  was  implicitly  applied  to 
the  residual  anomalies.  The  application  of  long-wavelength  terrain  correction 
to  approximate  a  first  order  long-wavelength  reduction  of  surface  anomaly  data 
to  a  level  surface  is  discussed  in  Moritz  (1966).  Since  the  final  position  of  the 
level  surface  to  which  the  data  are  reduced  is  uncertain  in  this  procedure,  it 
was  simply  assumed  that  this  position  coincides  with  the  mean  elevation  of 
topography  in  the  area  of  upward  continuation.  Such  uncertainties  in  defining 
the  data  level  directly  causes  uncertainties  in  defining  the  upward  continuation 
distance  H0  which  is  the  vertical  clearance  between  the  data  level  and  the 
upward  continuation  point  in  space.  The  numerical  effect  of  such  uncertainties 
on  upward  continuation  results  can  be  examined  from  Table  9. 

The  final  upward  continued  gravity  anomaly  model  of  the  indirect  method 
was  then  the  sum  total  of  three  contributions:  low  frequency  contribution  from 
spherical  harmonics,  high  frequency  contribution  from  topographic  mases,  and 
medium  frequency  contribution  from  Poisson  integration  of  terrain-corrected 
residual  field.  The  relative  order  of  magnitude  of  the  three  contributions 
depends  on  the  spectral  distribution  of  power  of  a  particular  field, and  for  the 
profiles  tested  in  Section  9  the  contributions  were  about  the  same  in 
magnitude. 

As  a  matter  of  interest,  we  compared  the  results  of  the  indirect  method  to 
those  of  two  simpler  upward  continuation  procedures.  The  first  was  the  direct 
Poisson  integration  of  the  original  terrain-uncorrected  surface  anomalies, 
simply  assuming  the  anomalies  to  lie  on  a  common  level  surface.  The  second 
was  the  direct  Poisson  integration  of  terrain-corrected  surface  anomalies,  with 
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the  terrain-corrections  intend  "d  to  effect  an  approximate  data  reduction  to  a 
level  surface.  The  terrain  corrections  in  the  second  method  originated  from 
the  point  terrain  corrections  given  in  the  original  data  records,  these 
corrections  having  been  applied  to  the  original  point  anomalies  before  the 
prediction  of  the  mean  anomalies  used  in  Poisson  integration. 

For  the  numerical  tests  we  developed  5’x5’  mean  anomalies  and  5'x5'  mean 
elevations  in  a  7*x9*  test  area  in  New  Mexico,  starting  from  available 
irregularly  distributed  point  anomaly  data  and  30"x30"  grid  point  elevations. 
The  30"x30"  elevations  themselves  were  also  used  in  the  final  computations,  for 
detailed  integration  of  topographic  effects  in  the  immediate  vicinity  of  the 
computation  point.  The  various  operational  steps  in  developing  mean  anomalies 
including  data  thinning,  tailoring  of  covariance  function,  use  of  only  the  ten 
closest  data  points  in  the  collocation  prediction,  and  data  de-trending  at  high 
and  low  frequencies  are  discussed  in  detail  in  Section  8.  The  required 
resolution  and  area  coverage  of  mean  anomalies  for  given  upward  continuation 
distances,  as  well  as  the  effect  of  data  error  propagation,  can  be  assessed 
based  on  concepts  presented  in  Sections  4,  5  and  6. 

Numerical  investigations  on  upward  continuations  to  te  >t-profiles  at  30,  10, 
and  5  km  are  presented  in  Section  9.  The  test  profiles  resulting  from  the 
direct  Poisson  integration  of  terrain-uncorrected  anomelies  are  negatively 
biased  (i.e.,  too  low)  by  about  (0.6,  0.5,  0.7)  mgal  at  elevation  (30,  10,  5)  km 
compared  with  the  profiles  resulting  from  the  direct  Poisson  integration  of 
terrain-corrected  anomalies.  This  bias  between  the  i  wo  direct  methods 
represents  the  effect  of  upward  continued  terrain  corrections  (see  (9.2.3)). 
There  is  no  detectable  bias  between  the  terrain-correctec  direct  method  and 
the  indirect  method;  this  is  mainly  due  to  the  fact  that  in  the  Poisson 
integration  part,  the  two  methods  both  use  terrain-corrected  anomalies  (see 
(9.2.6)  and  (9.2.7)).  The  standard  deviation  of  the  differences  among  all  three 
upward  continuation  methods  reach  the  order  of  (0.5,  0.6,  1.3)  mgal  at  (30,  10, 
5)  km  elevation  (see  Tables  3,  4,  and  5). 

In  Section  10  we  present  the  details  of  applying  the  upward  continuation 
methods  to  compute  anomalies  (and  total  gravity)  at  points  of  the  balloon-borne 
gravity  measuring  project  of  AFGL.  It  is  hoped  that  such  projects  would 
provide  "aerial  truth"  assessment  of  the  accuracies  of  upward  continuation 
models.  It  is  projected  in  Section  10.6  that  values  at  the  balloon  points  have 
been  recovered  with  about  0.9  mgal  standard  error  in  gravity  anomaly  with 
data  error  propagation  as  dominating  error  source,  and  about  0.7  mgal  error  in 
normal  gravity  with  vertical  position  error  as  dominating  error  source. 

In  another  series  of  tests  (Section  9.3)  we  have  shown  agreement  between 
Fast  Fourier  upward  continuation  and  Poisson  integration,  on  the  level  of  (0.1, 
0.3)  mgal  at  (30,  10)  km  elevation.  Fast  Fourier  techniques  are  useful  for  very 
fast  generation  of  complete  grids  of  upward  continued  values. 

We  conclude  that  actual  gravity  measurements  at  altitude  should  be  used  to 
validate  the  various  procedures  presented  in  this  report.  Practical  validation 
will  be  most  important.  First,  to  ascertain  the  accu  acy  of  the  upward 
continued  gravity  values,  and  second,  to  ascertain  the  improvements  gained  by 
employing  the  following  intended  refinements:  the  use  of  terrain  correction  to 
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DIGITAL  FOURIER  TRANSFORM  FOR  UPWARD  CONTINUATION 


//  JOB  , 

//  REGI0N=512K,MSGLEVEL={2»0) 

/* JOBPARM  V=S 

//PROCLIB  DD  D I SP=SHR»  DSN=GEODSC I .P ROCL I B 
//  EXEC  VSSUPER 
//SOURCE  DD  * 


UPWARD  CONTINUATION  OF  DG  BY  FOURIER  TRANSFORM 
IMPLICIT  REAL*8  (A-H.O-Z) 

C0MPLEX*16  Z( 128,128) ,CTEMP( 128) 

DIMENSION  TEMP (108) 

DIMENSION  01(2,128,128) 

EQUIVALENCE  (01(1,1,1)  ,ZI1,1)  ) 

COMMON  CTEMP 

INITIATE  CONSTANTS  (LINEAR  MEASURE  IN  KM  ANGULAR  INITIALLY  IN  DEC  DEG) 
H=28.5D0 
R =637 IDO *1.500 
AMELE V=1 • 5D0 
P I=4D0*0A  T  AN ( IDO) 

ANGRI 0=500/6000 
F I AVER=32. 500 

LINEAR  INCREMENTS  IN  E-W  AND  N-S  DIRECTIONS 

DX=R*OCCS(FIAVER*PI/180DC)*ANGRID*PI/18ODO 
OY=R*ANGRID*PI/ 18000 

INPUT  THE  DATA 
NX=128 
NY= 128 
INDX=108 
IN0Y=84 
C 

DO  2  1=1, NY 
DO  2  J= 1 , NX 
2  ZI I ,J)=( 000,000) 

DC  5  I = i » I NOY 

PEAD( 1,7001 )( TEMPI J) ,J=i ,INDX) 

DO  5  J= 1 ,  INOX 
Zl I ,J)=TEMP( J) 

5  CONTINUE 

7001  FORMAT ( 1  OX, l OF  7.2/1 3X,1 1F7.2)  ) 

C 


Y  _  r — y — j-  *  .  *.• — .▼T"  .  • 


C  FORWARD  2— C  TRANSFORM  (♦SIGN) 

CALL  FFT2D1Z,NY»NX,+1) 

C 

C 

C  FREQUENCY  INCREMENTS 
DFX=lDO/(NX*OX) 

DFY=100/(NY*DY) 

C 

C  UPWARD  CONTINUATION  2D  TRANSFORM  I  -  SIGN) 

CALL  FFT2UP1H,0FY,DFX,Z,NY,NX ,-l ) 

C - 

C 

C  OUTPUT 

DO  11  1*1 » INDY 
DO  10  J  =  l ,  I NDX 

10  TEMPI J)=Z1 I ,  J) 

11  IFd.EQ.42)  WRITE(6,7002)(  J, TEMPI  J)  ,J*l,INOX) 

7002  F0RMATI1X,I3»F7.2) 

C 

STOP 

END 

SUBROUTINE  FFT20  l  H, NX , NY, NS IGN > 

C  ♦♦  *********  *  *******  **********  ****.***  ****  ************* 

C  SUBROUTINE  FFT2D  COMPUTES  THE  TWO  DIMENSIONAL  FOURIER  TRANSFORM 

C  CF  A  COMPLEX  ARRAY  HINX,NY),NX  AND  NY  MUST  BE  A  POWER  OF  2. 

0 

C  NSIGN*  ♦!  INVERSE  TRANSFORM 

C 

C  NS  I GN=  -1  FORWARD  TRANSFORM 

C 

C* ********************************************* ******* 

IMPLICIT  REALMS  (A-H,0-Z) 

COMMON  CTEMP  . 

C0MPLEX*16  H(NX,NY),CTEMPU28) 

SIGNl=OFLOAT(NSIGN) 

00  10  IY=1 »NY 

10  CALL  FCRKINX,H(1,IY) .SIGN! ) 

IF(NY.EQ.l)  RETURN 
00  20  I X= 1 ,NX 
DO  30  IY*l,NY 
30  CTEMP I IY)=H( IX, IY) 

CALL  FCRK1NY, CTEMP, SIGNI  ) 

00  40  IV* l, NY 
40  H(IX, I Y)*CTEMP(IY) 

20  CONTINUE 
RETURN 
END 

SUBROUTINE  FOR K I  LX, CX , S IGN I ) 

£****« ************* ************* ********************** 

C  FAST  FOURIER  TRANSFORM,  MODIFIED  FROM  CLAERBOUT, J. F. 

C  FUNDAMENTALS  OF  GEOPHYSICAL  DATA  PROCESING,  MCCRAW-H I LL , l 976 . 

C  LX 

C  CX(K)=SUM(CX(J)*EXP(2*PI*SIGNI*I*(J-l)*(K-l)/LX)) 

C  J  =  1 
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C  FOR  K*l,2«***« ILX»2**INTEGER) 

C  S I GNI =  +1  INVERSE  TRANSFCRM 

C  SIGNI*  -I  FORWARD  TRANSFCRM 

C  LX  MUST  BE  A  POWER  OF  2  (LX»2** INTEGER  I 

C  NORMALIZATION  PERFORMED  EY  OIVIDING  BY 

C  DATA  LENGTH  UPON  THE  FORWARD  TRANSFORM 

£***** ******************************* ******************* 

IMPLICIT  REAL*8  (A-H,0-Zl 
COMPLEX+16  CX(LX) ,CARG,CEXP ,CW,CTEMP 
PI*4DO*DATAN( IDO) 

J=l 

SC" IOO/OFLOAT ( LX ) 

DO  30  1*1, LX 
IF(I.GE.J)  GOTO  10 
CTEMP*CXI J ) 

CXCJ)*CXI I  ) 

CXU)=CTEM  P 
10  M*LX/2 

20  IF1J.LE.M)  GOTO  30 
J*  J— M 
M=M/2 

IFIM.GE.I)  GOTO  20 
30  J= J+M 
L=  1 

40  I STEP=2*L 
OG  50  M*  1 ,  L 

CARG*I0D0,1D0)*(  P I*S  IGNI  ♦DFLOAT  (M-l  D/DFLOAT  (L  ) 

CW*CDE  XP ( C ARG ) 

DO  50  I*M,LX, I  STEP 
CTEMP=CW*CX 1 1 >L ) 

CXI I+L)*CX( IJ-CTEMP 
50  CX(I)*CXm+CTEMP 
L* I  STEP 

IF(L.LT.LX)  GOTO  40 
IFISIGNI.GT.ODOI  RETURN 
00  60  1*1, LX 

60  cxm=cxm*sc 

RETURN 
ENO 

SUBROUTINE  FFT2UP  ( ELEV  ,  CFX , DFY , H, NX , NY , NS IGN ) 
***************************************************** 

SUBROUTINE  FFT2D  COMPUTES  THE  TWO  DIMENSIONAL  FOURIER  TRANSFORM 
CF  A  COMPLEX  ARRAY  H(NX,NY),NX  AND  NY  MUST  BE  A  POWER  OF  2. 
ELEV  IS  THE  APWARD  CONTINUATION  DISTANCE  l SAME  UNITS  AS  DX  DY) 

DFX  DFY  FREQUENCY  INCREMENTS  ASSOCIATED  WITH  INDEXES  NX  NY  RESP 
DFX=l/INX*DX)  DFY*1/ I NY*DY ) 

NS  IGN*  *1  INVERSE  TRANSFCRM 

NSIGN*  -1  FCRWARO  TRANSFCRM 

***************************************************** 
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IMPLICIT  REAL* 8  (A-H.C-Z) 

COMMON  CTEMP 

C0MPLEX*16  H<NX,NY) ,CTEMP( 128) 

PI=4DO*OATAN( 100) 

S16NI=DFL0AT(NSIGN) 

C 

DO  16  I  X=  1 , NX 

IF ( IX- l.GT  .NX/2 ) THEN 

IPX=I IX— 1 )-NX 

ELSE 

IPX=IX-1 

END  IF 

DO  15  IY-l.NY 

IF1IY-1.GT.NY/2)  THEN 

IPY=I  IY-D-NY 

ELSE 

IP Y=I Y- 1 

ENOIF 

DUMPED EX P { -ELEV*2D0*P I *D  SORT ( DFLO AT (IPX )  **2*0FX**2 ♦□FLOAT < I P Y ) 
$**2*DFY**2) ) 

15  H( IX, IY)=DUMP*H( IX, IY) 

16  CONTINUE 
C 

DC  10  I Y  =  1 » NY 

10  CALL  FCRK(NX,H(1,IY)  ,SIGM  ) 

IF(NY.EQ.l)  RETURN 
C 

DO  20  I X®1 , NX 
DO  30  I Y—  1 ,  NY 
30  CTEMPt  IY)«M  IX.IY) 

CALL  FORK ( NY, CTEMP , S IGNI  ) 

DO  40  I Y= 1 , NY 
40  H(IX,IY)=CTEMP(  IY) 

20  CONTINUE 
RETURN 
END 

/* 

//GO.FTOIFOO  1  DO  D I  SP-=SHR ,  0SN=TS426 8. FA5X 58 .NEWMEX 
// 
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